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ABSTRACT 

A long-standing problem in low-mass star formation is the "luminosity problem," whereby proto- 
stars are underluminous compared to the accretion luminosity expected both from theoretical collapse 
calculations and arguments based on the minimum accretion rate necessary to form a star within the 
embedded phase duration. Motivated by this luminosity problem, we present a set of evolutionary 
models describing the collapse of low-mass, dense cores into protostars. We use as our starting point 
the evolutionary model following the inside-out collapse of a singular isothermal sphere as presented by 
Young & Evans (2005). We calculate the radiative transfer of the collapsing core throughout the full 
duration of the collapse in two dimensions. From the resulting spectral energy distributions, we cal- 
culate standard observational signatures (Lf,oi, Ti,oi, L^oi / L smm) to directly compare to observations. 
We incorporate several modifications and additions to the original Young & Evans model in an effort 
to better match observations with model predictions: (1) we include the opacity from scattering in the 
radiative transfer, (2) we include a circumstcUar disk directly in the two-dimensional radiative trans- 
fer, (3) we include a two-dimensional envelope structure, taking into account the effects of rotation, 
(4) we include mass-loss and the opening of outflow cavities, and (5) we include a simple treatment of 
episodic mass accretion. We find that scattering, two-dimensional geometry, mass-loss, and outflow 
cavities all affect the model predictions, as expected, but none resolve the luminosity problem. On 
the other hand, we find that a cycle of episodic mass accretion similar to that predicted by recent 
theoretical work can resolve this problem and bring the model predictions into better agreement with 
observations. Standard assumptions about the interplay between mass accretion and mass loss in our 
model give star formation efficiencies consistent with recent observations that compare the core mass 
function (CMF) and stellar initial mass function (IMF). Finally, the combination of outflow cavities 
and episodic mass accretion reduce the connection between observational Class and physical Stage to 
the point where neither of the two commonly used observational signatures (Tboj and Lhoi/ Lsmm) can 
be considered reliable indicators of physical Stage. 
Subject headings: stars; formation - stars: low-mass, brown dwarfs 



1. INTRODUCTION 

Over the past few decades a general picture of low- 
mass star formation has emerged. As first presented by 
Adams et al. (1987) and summarized by Shu, Adams, & 
Lizano (1987), this picture merges an empirical classifi- 
cation scheme based on the infrared spectral slope (Lada 
& Wilking 1984) with a theory involving the stages of 
the collapse of a dense, rotating core (Shu 1977; Tere- 
bey, Shu, & Cassen 1984, hereafter TSC84). In Stage 
I, the core begins collapsing and the newly formed pro- 
tostar^ is initially heavily obscured by the surrounding 
envelope, exhibiting a Class I spectral energy distribu- 
tion (SED) rising from 2 to 20 fim due to reprocessing 
of short-wavelength emission by the dust in the enve- 
lope. Conservation of angular momentum causes a disk 
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to build up (e.g., Adams & Shu 1986). The envelope dis- 
sipates through accretion and mass-loss processes. Once 
it fully dissipates the object transitions from Stage I to 
Stage II, leaving a pre-main sequence star surrounded by 
a circumstellar disk that exhibits a Class II SED falling 
from 2 to 20 /xm, but with a shallower slope than ex- 
pected for a main-sequence star due to "extra" infrared 
emission from the dust in the disk. The disk eventu- 
ally dissipates, leaving a Stage III pre-main sequence 
star exhibiting a Class III SED consistent (or at least 
nearly so; see Evans et al. [2009] and references therein) 
with that expected for a main-sequence star. Andre et 
al. (1993) later added Class to the scheme, defining such 
objects observationally as emitting a relatively large frac- 
tion (greater than 0.5%) of their total luminosity at wave- 
lengths A > 350 /im. Defining a corresponding physical 
stage. Stage objects are young, embedded protostars 
with greater than 50% of their total system mass still in 
the envelope (Andre et al. 1993). While the term "Class" 
is often assumed to have both meanings in the literature, 
in this work we follow Robitaille et al. (2006) and distin- 
guish between "Class" , determined by observed quanti- 
ties, and "Stage" , determined by the ratio of envelope 
mass to total system mass. 

Despite the successes of this picture many questions 
remain, including a detailed understanding of the mass 
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accretion process from the core to the star. The "stan- 
dard model" of star formation, the inside-out collapse 
of an isothermal sphere calculated by Shu (1977) and ex- 
tended by TSC84 to include rotation, predicts a constant 
mass accretion rate of about 2 x 10"^ M0 yr~^. This 
gives rise to the classic "luminosity problem" whereby 
accretion at such a rate produces accretion luminosities 
{Lace oc M^M) higher than typically observed for embed- 
ded protostars. First noticed by Kenyon et al. (1990), 
the problem has recently been emphasized by results 
from the Spitzer Space Telescope. Dunham et al. (2008), 
Enoch et al. (2009), and Evans et al. (2009) all show 
that the distribution of embedded protostellar luminosi- 
ties is strongly peaked at low luminosities. Enoch et 
al. and Evans et al. both find that a substantial fraction 
(greater than 50%) of embedded protstars have luminosi- 
ties suggesting accretion rates M < 10~^ Mq yr~^, and 
Dunham et al. argue that the large fraction of sources 
at low luminosities is inconsistent with a constant mass 
accretion rate. 

To compare theoretical models of star formation to ob- 
servations. Young & Evans (2005; hereafter YE05) used 
a one-dimensional dust radiative transfer package to cal- 
culate the observational signatures of cores undergoing 
inside-out collapse following Shu (1977). They followed 
three different cores with initial masses of 0.3, 1, and 3 
Mq from the onset of collapse until the end of the embed- 
ded phase, calculating the bolometric luminosity (Ltoi), 
bolometric temperature {Thoi), and ratio of bolometric to 
submillimeter luminosity {Lhoi/ Lsmm, see ij3.1[) . Thoi is 
defined by Myers & Ladd (1993) as the temperature of a 
blackbody with the same flux-weighted mean frequency 
as the source (see §3.1|1 . and can be thought of as a pro- 
tostellar equivalent of T^f f ; Tboi starts at very low values 
10 K) for cold, starless cores and eventually increases 
to Teff once the envelope and disk have fully dissipated. 
YE05 compared their model to observations by plotting 
both their model tracks and observations of sources on 
a plot of Lboi vs. Tfjoi, which Myers et al. (1998) called 
a BLT diagram. This figure (Figure 19 in YE05) shows 
that observed sources at a given Tboi range from having 
Lboi consistent with the Young & Evans model tracks to 
having Lhoi up to 1 — 3 orders of magnitude lower than 
predicted, clearly illustrating the luminosity problem. 

An idea proposed to resolve the luminosity problem 
is that mass accretion is episodic in nature, and the 
protostars with the lowest luminosities are those ob- 
served in quiescent accretion states (e.g., Kenyon et 
al. 1990, Kenyon & Hartmann 1995; YE05; Enoch et 
al. 2009; Evans et al. 2009). Theoretical studies have 
provided several mechanisms by which such a process 
may occur, such as material piling up in a circumstel- 
lar disk until gravitational instabilities drive angular mo- 
mentum outward and mass inward in short-lived bursts 
(Vorobyov & Basu 2005b, 2006; Boss 2002). Ahcrna- 
tively, accretion bursts may be driven by a combination 
of gravitational and magnetorotational instabilities (Zhu 
et al. 2009), or quasi-periodic magnetically driven out- 
flows in the envelope may cause mass accretion onto the 
protostar to occur in magnetically controlled bursts (Tas- 
sis & Mouschovias 2005). Indeed, observational evidence 
for non-steady mass accretion in young protostellar sys- 
tems still in the embedded phase now exists in the form of 



accretion bursts in Class I sources (e.g., Acosta-Pulido et 
al. 2007; Kospal et al. 2007; Fedele et al. 2007) and Class 
sources with strong outflows implying higher average 
mass accretion rates than expected from currently ob- 
served low luminosities (e.g., Dunham et al. 2006; Andre 
et al. 1999; M. M. Dunham et al. 2009, in preparation). 
Additionally, Watson et al. (2007) showed a mismatch 
between the accretion rates onto the disk and protostar of 
NGC 1333-IRAS 4B (measured by modeling water emis- 
sion lines and by assuming all of the observed luminos- 
ity is accretion luminosity, respectively), a result they 
have now expanded to other sources (D.M. Watson et 
al. [2009], in preparation). Finally, episodic mass ejec- 
tion is seen in jets ejected from some protostellar sys- 
tems, suggesting an underlying variability in the mass 
accretion, although the combination of jet velocities and 
spacing between knots often suggests shorter periods of 
episodicity than found by the above theoretical studies. 
For example, Lee et al. (2007) found a period of ~ 15 — 44 
yr for the periodic protostellar jet HH 211. We also note 
here that an alternative collapse scenario, "outside-in" 
collapse, where the collapse is triggered by an external 
shock wave, can produce a range of mass accretion rates 
roughly in agreement with those derived from observa- 
tions and predicted by episodic accretion models (Boss 
1995). However, while such a process is relevant for mas- 
sive star-forming regions and possibly for the formation 
of our own solar system (Boss 2008), it is not likely rele- 
vant for the relatively isolated protostars forming in most 
nearly, low-mass star forming regions. 

Here we will test the hypothesis that episodic accre- 
tion can solve the luminosity problem. First, however, 
we will test the effects of several other possibilities that 
were not included in the YE05 model, including revised 
dust opacities, two-dimensional disk and envelope geom- 
etry, and mass- loss and outflow cavities. This work is 
complementary to several other recent modeling efforts. 
Myers et al. (1998) included the effects of mass-loss in 
their calculations of the evolution of Lf,oi and Ttoi with 
time for collapsing cores, but they did not include out- 
flow cavities and their model evolution is not based on 
a fully self-consistent model such as the collapse solu- 
tions calculated by Shu (1977) or TSC84. Whitney et 
al. (2003a, 2003b), Robitaille et al. (2006), and Crapsi 
et al. (2008) all used 2-D radiative transfer models to 
investigate the effects of two-dimensional disk and enve- 
lope geometry and outflow cavities on the evolutionary 
signatures of embedded protostars. However, none of 
these authors present a full evolutionary model following 
the collapse of a core but instead vary parameters over 
pre-defined grids to capture typical protostars of differ- 
ent evolutionary stages, and only Crapsi et al. (2008) 
considered the predictions of observed quantities other 
than infrared colors. Lee (2007) included episodic accre- 
tion in the YE05 model in a very simple manner in order 
to study the effects such a process has on the chemical 
evolution of collapsing cores. Myers (2008) presented an 
analytic calculation of the growth of a protostar through 
competing infall and dispersal processes; some aspects of 
our models, in particular the opening of outflow cavities, 
are similar to those featured by Myers (2008). Baraffe, 
Chabrier, & Gallardo (2009) showed that episodic ac- 
cretion in the early, embedded phase can explain the 
observed luminosity spread in H— R diagrams of star- 
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forming regions at a few Myr without having to invoke 
large age spreads. FinaUy, Vorobyov (2009b) compared 
the distribution of mass accretion rates in their simu- 
lations of episodic accretion (Vorobyov & Basu 2005b, 
2006) to those inferred from the luminosities of proto- 
stars in the Perseus, Serpens, and Ophiuchus molecular 
clouds compiled by Enoch et al. (2009) and showed that 
their simulations reproduced some of the basic features 
of the observed distribution of mass accretion rates. 

In this paper we revisit the YE05 model, which is sum- 
marized in 321 Following YE05, we assume a distance of 
140 pc for all models to calculate observed SEDs. This 
assumed distance only affects the absolute flux level when 
we display SEDs, all other observational signatures are 
independent of the assumed distance. We discuss the 
method we use to compare evolutionary models to ob- 
servations in Sj3l In 21 we make several updates and ad- 
ditions to the model in a step-by-step fashion, examining 
the results of each modification individually. Specifically, 
we include scattering in the radiative transfer calcula- 
tions in In i j4. 21 - 14.31 we generalize the model from 
its original, one-dimensional structure to two dimensions 
with a more realistic disk (i i4.2|) and two-dimensional en- 
velope structure ( §4.3|) . We include the effec ts of mass- 
loss and outflow cavities in m.4[ and in H4.5I we include 
a simple treatment of episodic accretion. A discussion of 
the results of this work is presented in ^ and we present 
a summary of our conclusions in Sj6l We note here that 
choices of parameters in the models presented below are 
physically motivated and theoretically and/or observa- 
tionally constrained whenever possible. However, these 
are simple, idealized models that are not always fully 
self-consistent. Limitations are discussed as each modi- 
fication is described. 

2. DESCRIPTION OF THE ORIGINAL MODEL 

We begin with a summary of the YE05 model. We 
provide a fairly comprehensive description of this model 
to place our work in later sections in context, but refer 
the reader to YE05 for a complete description. 

2.1. Evolution of the Envelope, Protostar, and Disk 

The YE05 model follows the collapse of singular 
isothermal spheres with initial masses of 0.3, 1, and 3 
M0 according to the inside-out collapse solution calcu- 
lated by Shu (1977). This model begins with an envelope 
radial density profile proportional to r'^^^, truncated at 
an outer radius that sets the initial core mass (YE05 
Equation 1)''. The collapse of the envelope begins at the 
center and moves outward at the sound speed, giving 
rise to an infall radius that moves outward with time. 
The density profile is then described approximately as 

3/2 

a broken power law; inside the infall radius n oc renv , 
indicative of free-fall, while outside the infall radius the 
profile remains the initial n cx r~^^ . YE05 used the exact 
solutions from Shu (1977) to account for the transition 
region between the two. Once the infall radius exceeds 
the envelope outer radius, the model adopts a density 

3/2 

profile with n cx Venv everywhere. The envelope inner 

Following the convention adopted by YE05, radii pertaining to 
the envelope are denoted by lowercase r, while radii pertaining to 
either the star or disk are denoted by uppercase R. 



radius is held fixed at a value such that the initial optical 
depth at 100 /xm is set equal to 10 (YE05 Equation 4; 
see YE05 for a discussion of the effects of varying this 
initial optical depth) until the disk outer radius exceeds 
this value (see below) ; once this occurs the inner envelope 
radius is set equal to the disk outer radius. 

The mass of the envelope, M^nv , declines as mass ac- 
cretes from the envelope to the protostar -|-disk system at 
the rate Menv In the Shu (1977) solution, this rate is 
constant and given by 



M, 



env — ^0 Q 



where ttiq is a dimensionless constant of order unity, G is 
the gravitational constant, and Cs is the effective sound 
speed including thermal and turbulent components and 
calculated by YE05 as 
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With k the Boltzmann constant, T the isothermal tem- 
perature (assumed to be 10 K), /i the mean molecular 
mass (2.29 for a gas that is 25% by mass helium), niH 
the mass of the hydrogen atom, and Vturb the turbulent 
velocity, chosen such that the thermal and turbulent con- 
tributions are equal, Cg = 0.268 km s~^. This gives an 
envelope mass accretion rate of 4.57 x 10~^ Mq yr~^. 
This is 5% lower than the YE05 value of Menv = 4.8 x 
10~^ Mq yr~^; this difference arises from correcting a 
small numerical error in the YE05 model. 

To include a protostellar disk in their 1-D model, YE05 
use the method developed by Butner et al. (1994), based 
on the disk model of Adams et al. (1988). This method 
simulates a disk by calculating the emission from a disk 
with given surface density and temperature profiles at 
a given inclination, averaging over all inclinations, and 
then adding this average spectrum to the (proto) stellar 
spectrum to form the final input spectrum of the internal 
source for the 1-D radiative transfer calculation through 
the envelope. Both the surface density and temperature 
profiles are described as power laws: T,{Rdisk) oc -R^j^j, 
and T{Rdisk) oc R^ig^.- YE05 choose p = 1.5, following 
Butner et al. (1994) and Chiang & Goldreich (1997), and 
q = 0.5 to simulate a flared disk. 

The inner radius of the disk is set equal to the dust de- 
struction radius, calculated (assuming spherical, black- 
body dust grains) as 



R 



disk 
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where is the protostellar luminosity (see below) and 
Tdust is the dust destruction temperature. YE05 assume 
Tdust = 2000 K; here we update this value to Tdust = 
1500 K (e.g., Cieza et al. 2005). The outer radius of the 
disk is set to the centrifugal radius, which increases with 
time as (TSC84) 



^disk 
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(4) 



where t is the time since the onset of the collapse and 51o 
is the initial angular velocity of the cloud. YE05 set fio 
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such that the disk outer radius is 100 AU at the end of 
the collapse of each core (fla = 3.4 x 10"^^, 5.5 x 10~^^, 
and 1.0 X 10"" s'^ for the 0.3, 1, and 3 M© cores, re- 
spectively). The 5% lower envelope accretion rate from 
YE05 results in 5% longer total collapse duration, but 
since wc choose to keep the values of fio set by YE05 
to minimize changes and facilitate direct comparison be- 
tween their results and the results of this work, the final 
disk outer radii are approximately 15% larger than 100 
AU. 

Following Adams fc Shu (1986), YE05 assume that 
all mass accreted from the envelope accretes onto ci- 
ther the star or the disk at rates M* and Mdisk, where 
Menv = Af* + Mdisk- These rates are governed by 
M*, the ratio of the protostellar and disk outer radii 
(u* = i?*/i?2"jj,). With the protostellar radius calcu- 
lated following Palla & Stabler (1991), except at early 
times {t <2x 10^ yr) where it is set to 5 AU to simulate 
the First Hydrostatic Core (Masunaga et al. 1998; Boss 
& Yorkc 1995; sec YE05 for details), Adams & Shu use 
the velocity field and density profile of the collapse so- 
lution for a rotating, singular isothermal sphere (Cassen 
& Moosman 1981; TSC84) to determine the flux of ma- 
terial flowing from the cloud directly onto the protostar 
and disk and calculate the protostellar and disk mass 
accretion rates as 

=Me„.[l-(l-w*)'/'] , (5) 

Md^sk ^ Menvil ~ U,y/'^ . (6) 

In addition to direct accretion from the envelope 
(which quickly becomes negligible as the disk outer ra- 
dius grows and thus decreases), the protostar also ac- 
cretes mass from the disk at the rate MjjtoP = VoMenv, 
where rjjj is an efficiency factor assumed to be 0.75 (see 
YE05 for a discussion of the effects of different assumed 
values for ??£)). Following Adams & Shu (1986), YE05 
calculate the mass of the disk, Mdisk, including both ac- 
cretion from the envelope and onto the protostar (see 
YE05 Equations 12 and 13). The mass of the proto- 
star is then calculated as = Mj„t — Mdisk, where 
Mint is the total internal mass accreted from the enve- 
lope {Mint = Menvt)- 

2.2. Luminosity Sources 

The total internal luminosity of the protostar and disk 
at each point in the collapse from core to star contains 
several components. Following Adams & Shu (1986), 
YE05 include six components: 

1. LEtoP- luminosity arising from accretion from the 
envelope directly onto the protostar. 

2. LEtoD- luminosity arising from accretion from the 
envelope onto the disk. 

3. LntoP'- luminosity arising from accretion from the 

disk onto the protostar. 

4. LjjM'- disk "mixing luminosity" arising from lu- 
minosity released when newly accreted material 
with its own radial and angular velocity compo- 
nents mixes with disk material in a Keplerian orbit. 



putting the new and old material into a new Kep- 
lerian orbit (see Adams & Shu [1986] for details). 

5. Ljjn: luminosity arising from the release of energy 
stored in differential rotation of the protostar. 

6. Lphot- luminosity arising from gravitational con- 
traction and deuterium burning. 

The first two components are both directly propor- 
tional to Menv, with geometrical factors that depend on 
to account for the changing rates of accretion onto the 
protostar and disk (see Equations 27 and 29a of Adams 
& Shu 1986 for the exact definitions of each term). Both 
components are quite small throughout the collapse of 
each core; LeloP since the amount of material accret- 
ing directly from the envelope to the star becomes small 
very quickly, and LEtoD because -R^isfe >> R* (except 
for very early times) and thus the material accreting onto 
the disk has not yet fallen very deep into the potential 
well. 

The third component depends on the rate at which ma- 
terial accretes from the disk to the protostar, controlled 
by the efficiency factor rjo- The exact definition is given 
by Equation 30 of Adams & Shu (1986); it is essentially 
just one-half of the spherical accretion luminosity arising 
from material accreting at this rate (the other half of the 
initial gravitational potential energy is stored as kinetic 
energy of the disk material's Keplerian rotation shortly 
before it accretes onto the star), with both the luminosity 
already released from accretion onto the disk and from 
the disk mixing (see below) subtracted out. This is the 
dominant source of luminosity once the disk has formed. 
Any luminosity arising from the inward transport of ma- 
terial within the disk is indirectly included in this term 
since its calculation starts with the total spherical accre- 
tion luminosity. 

The fourth component is also directly proportional to 
Menv, with geometrical factors that depend on u*. The 
exact definition is given by Equation 29b of Adams & 
Shu (1986). 

The fifth component depends on the total rate of ac- 
cretion onto the star and the efficiency ?7# with which 
energy stored in differential rotation of the protostar is 
released (assumed to be 77, = 0.5; see YE05 for details). 
The exact definition is given bv Equation 32 of Adams 
& Shu (1986). 

To include the sixth component, the luminosity arising 
from gravitational contraction and deuterium burning, 
YE05 used the pre-main sequence tracks of D'Antona 
& Mazzitelli (1994) with opacities from Alexander et 
al. (1989). They also assumed a power-law expression 
to extrapolate to times earlier than those included in the 
pre-main sequence tracks, Lphot = L^°*{t/to)^ , where to 
is the earliest time in the tracks and L^°^ is the pre-main 
sequence luminosity at this time. Finally, they followed 
Myers et al. (1998) in adding 10^ yr to the times of the 
pre-main sequence tracks to account for the delay be- 
tween the onset of collapse and the "zero time" of these 
tracks. Thus it is only at late times in the collapse of the 
cores that Lphot becomes an important source of lumi- 
nosity. 

Finally, there is also external luminosity arising from 
heating of the envelope by the Interstellar Radiation 
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Field (ISRF). We adopt the same ISRF as YE05, and 
input the mean intensity of the ISRF (Jext) into the ra- 
diative transfer code as an additional source of heating. 
The luminosity added to Lboi from this external heating, 
Lg^t, is determined after each radiative transfer model 
is run by subtracting the total internal luminosity (the 
sum of the above six components) from the total model 
luminosity. 

3. COMPARING MODELS TO OBSERVATIONS 

For all of the models presented below, we use the 
two-dimensional, axisymmetric, Monte Carlo dust ra- 
diative transfer package RADMC (DuUemond & Turolla 
2000; DuUemond & Dominik 2004) to calculate the two- 
dimensional dust temperature profile of the envelope at 
each timestep throughout the collapse of the 0.3, 1, and 
3 Mq cores (following YE05, At = 1000, 2000, and 6000 
yr for the 0.3, 1, and 3 Mq cores, respectively). Spec- 
tral Energy Distributions (SEDs) at each timestep are 
then calculated at 9 different inclinations ranging from 
i = 5 - 85° in steps of 10° except for model 1 (SU, 
where there is no inclination dependence and thus only 
one SED is calculated per timestep. An inclination of 
i = 0° corresponds to a pole-on (face-on) system, while 
an inclination of i = 90° corresponds to an edge-on sys- 
tem. 

3.1. Calculating Observational Signatures 

We use the model SEDs to calculate observational 
signatures of the models at each timestep for each 
inclination. We calculate the bolometric luminosity 
(Lboi), the ratio of bolometric to submillimeter luminos- 
ity {Lboi / L smm) , and the bolometric temperature {Tioi)- 
Lboi is calculated by intergrating over the full SED, 

/•OC 

Lbol = ind^ / S^di^ , (7) 
Jo 

while the submillimeter luminosity is calculated by inte- 
grating over the SED for A > 350 ^m, 

/"Z^— c/350 /^m 

Lsmm = 47rd^ / Si.dv . (8) 

Jo 

The bolometric temperature is defined to be the temper- 
ature of a blackbody with the same flux-weighted mean 
frequency as the source (Myers & Ladd 1993). Following 
Myers & Ladd, Tboi is calculated as 

Tboi = 1.25 X 10-" -^"oo ^ , K . (9) 
Jo ^'^'^^ 

The integrals defined in equations [7| — [9] are calculated 
using the trapezoid rule to integrate the finitely sampled 
model SEDs. 

Both Tfjoi and Lboi/ Lsmm can be used as alternatives 
to the infrared spectral slope to classify sources. Evans 
et al. (2009) present a comprehensive discussion of source 
classification by observational signatures, here we briefly 
summarize the main points. Chen et al. (1995) defined 
class boundaries in Tboi, as follows: 

Class Tboi < 70 K 

Class I 70 < Tboi < 650 K 



Class II 650 < Tboi < 2800 K 

From the original observational definition of Class by 
Andre et al. (1993), the class boundaries in Lboi/Lsmm 
are: 

Class Lbol /Lsmm < 200 

Class I Lbol /Lsmm > 200 

Based on their evolutionary model, YE05 revised the 
Class O/I boundary in Lboi/Lsmm from 200 to 175. There 
is no defined boundary between Class I and Class II in 

Lbol / Lsmm- 

3.2. Observational Dataset 

We use the 1024 Young Stellar Objects (YSOs) in 
the five large, nearby molecular clouds surveyed by the 
Spitzer Space Telescope Legacy Project "From Molecular 
Cores to Planet Forming Disks" (Evans et al. 2003) as our 
observational dataset. Evans et al. (2009) compiled pho- 
tometry and calculated Lboi and Tboi for all 1024 YSOs 
in the same manner as described above. They used both 
the observed photometry to calculated observed Lboi and 
Tboi values, and photometry corrected for foreground ex- 
tinction to calculate extinction-corrected values of Lboi 
and Tboi (see Evans et al. [2009] , and ^5.'2[ for details on 
the corrections for foreground extinction). Since evolu- 
tionary models have no foreground extinction, only local 
extinction from the dust in the model itself, we use the 
extinction-corrected Lboi and Tboi. 

Evans et al. (2009) concluded that 112 of the 1024 
YSOs are embedded protostars based on association with 
a millimeter continuum emission source tracing an enve- 
lope. Since the evolutionary models presented both by 
YE05 and in this paper follow the evolution of protostars 
up until the point of complete envelope dissipation, we 
consider these 112 embedded protostars to be our final 
observational dataset to use when comparing the models 
to observations. 

3.3. Comparing Models to Observations 

With the model observational signatures calculated as 
described above and the observational dataset of 112 em- 
bedded protostars from Evans et al. (2009), we can com- 
pare the model predictions to observations. 

The most common method of comparing model pre- 
dictions to observations is by plotting the model tracks 
on a diagram of Lboi vs. Tboi- This was first done by 
Myers et al. (1998), who called such a diagram a BLT 
diagram. Figure [1] shows the YE05 model tracks for the 
0.3, 1, and 3 Mq cores plotted on a BLT diagram, similar 
to Figure 19 from YE05. Also plotted arc the 1024 YSOs 
from Evans et al. (2009), with color indicating spectral 
class (red for Class O/I, green for Flat spectrum, blue for 
Class II, and purple for Class III; see Evans et al. for 
details) and symbol indicating source type (circles for 
sources associated with envelopes as traced by millimeter 
continuum emission, plus signs for sources not associated 
with envelopes) . The relevant comparison is between the 
model tracks and the sources plotted as circles. 

Plotting Lbol — Tboi tracks on a BLT diagram such as 
in Figure [1] is an adequate means of comparing model 
results to observations for the YE05 model. There is 
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Fig. 1.— BLT diagram for the YE05 model. The Ught blue 
lines show Li,oj — Tj,oj tracks for the 0.3, 1, and 3 Mq cores. The 
grayscale pixels indicate the fraction of total tim e the model spends 
in each L^oj — T;,;,; bin, calculated from Equation llOl The grayscale 
is displayed in a logarithmic stretch with the scaling chosen to 
emphasize the full extent of the models in L^o; — Ti,ai space. The 
mapping between grayscale and fraction of total time is indicated 
in the legend. The class boundaries in T^o; are taken from Chen 
et al. (1995). The thick dashed Une shows the ZAMS (D'Antona 
& MazziteUi 1994) from 0.1 to 2.0 Mq. The colored symbols show 
the Young Stellar Objects from Evans et al. (2009) in this diagram; 
the color indicates spectral class (red for Class O/I, green for fiat 
spectrum, blue for Class II, and purple for Class III), while a circle 
or cross indicates the source is or is not, respectively, associated 
with an envelope as traced by millimeter continuum emission. 

no inclination dependence so there is only one track per 
core mass, and both Lboi and Tboi increase monotonically. 
However, most models below will have an inclination de- 
pendence, making it difficult to compare model tracks 
to observations when each model has different tracks de- 
pending on inclination. Furthermore, the change in L^oi 
and Thoi with time will no longer be monotonic once 
episodic accretion is introduced ( ii4.5[) . eliminating the 
concept of "model track" altogether. 

With this in mind, we introduce other methods of com- 
paring to observations. First, we divide the Ltoi — T^oi 
space into bins of 1/3 dex in both dimensions. We then 
calculate the fraction of total time the model spends in 
each Lboi - T^oi bin {fun) as follows: 




where the numerator is the time spent in the bin and the 
denominator is the total time. The interior sum in both 
the numerator and denominator is over the 9 different 
inclinations while the exterior sum is over the three dif- 
ferent initial core masses. The quantity tine is the total 
time that the SED at a given inclination spends in the 
specified L^oi ~ T^oi bin whereas tcoUapse is the total col- 
lapse time of the core (67,000, 224,000, and 673,000 yr 
for the 0.3, 1, and 3 M© cores, respectively). Winc is the 
weight each inclination receives in the sum, defined as 
the fraction of solid area subtended by that inclination. 
This is calculated in practice by assuming each of the 
9 SEDs calculated is valid for inclinations spanning the 
range (i - 5°) to (i-l-5°). 

The final quantity in Equation [10] is Wmass, the weight 
given to each of the three initial mass cores. This is de- 



termined by the core mass function (CMF) of starless 
cores. Reading directly from the CMF plot presented 
by Enoch et al. (2008; their Figure 13), we find a ratio 
of 1 to 3 and 1 to 0.3 Mq cores of N1/N3 = 2.3 and 
Ni/Nq,3 = 2.8. Alternatively, using their best-fit power- 
law of dN/DM cx M-2-3±o-4 givesS N1/N3 = 3.91:1a, 
while using their best-fit lognormal distribution gives 
N1/N3 = 1.1 and Ni/No.3 = 12.9. If we instead use 
the three-component power-law IMF found by Kroupa 
(2002) and assume a 30% star-formation efficiency in 
dense gas (Alves et al. 2007; see also i 34.4p to scale 
from the IMF to the CMF, we obtain N1/N3 = 2.3 and 
Ni/No,3 = 0.6. Finally, if we assume the same efficiency 
but instead use the IMF found by Muench et al. (2002) 
for the Trapezium cluster shown by Alves et al. (2007) 
to match (with the 30% scaling) the dense core mass 
function in the pipe nebula, we find N1/N3 = 1.1 and 
Ni/No,3 — 3.5. Given the uncertainties in determining 
the exact CMF, both from uncertainties in core mass 
calculations and from completeness effects (see Enoch 
et al. 2008 for a complete discussion), we simply av- 
erage the above values^ to obtain N1/N3 = 2.1 and 
Ni/Nq,3 = 2.3. Requiring the sum of the weights to 
be 1 gives wo.3 = 0.2275, wi = 0.5233, and W3 = 0.2492. 

In addition to the model tracks. Figure [T] also shows, 
using grayscale pixels, the fraction of total time the YE05 
model spends in each Lhoi — Tboi bin, calculated from 
Equation [10] above. Comparison to the model tracks il- 
lustrates that this method of displaying the results has 
the advantage of showing not only what regions of the 
BLT diagram the model encompasses but also the rel- 
ative amount of time spent in different portions of the 
diagram. For all of the revised models presented in this 
paper, we do not show model tracks and instead use only 
this method of displaying the results on the BLT dia- 
gram. 

We can also compare the overall distribution of time 
the models spend at different values of Ltoi and Ttoi to 
the fraction of total sources observed at those values. As 
an example. Figure [2] shows Li,oi and Ttoi histograms for 
both the observations and the YE05 model, with bin- 
sizes of 1/3 dex for both quantities. The observational 
histograms only include the 112 embedded sources as- 
sociated with envelopes (plotted as filled circles on the 
BLT diagrams) and plot the fraction of total sources in 
each bin, while the model histograms plot the fraction 
of total time spent in each bin calculated from Equation 
[TUl This figure emphasizes the higher luminosities of the 
model compared to the observations. We can quantify 
the agreement between the model and the observations 
with K-S tests. Such tests shows that there is less than 
a 0.1% probability that the observed and model Lboi his- 
tograms represent the same underlying distribution, and 
a 56% probability that the observed and model Tboi his- 
tograms represent the same underlying distribution. 

Finally, we also divide the BLT diagram into much 
larger bins (1 dex in Tboi and 2 dex in Lboi', the bins are 

* As the power-law is only fit to the CMF for M > 0.8 Mq, it 
is inappropriate to use it to obtain an estimate of N1/N0.3 

^ We leave Ni/Nq,3 = 12.9, obtained from the best-fit lognormal 
distribution in Enoch et al. (2008), out of the average since it is sig- 
nificantly higher than other values and is derived from a lognormal 
fit to data that likely suffers from incompleteness effects. 
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Fig. 2. — Histograms showing the fraction of total sources (ob- 
servations; soUd lines) and fraction of total tim e sp ent by the YE05 
model (dashed lines; calculated from Equation 110^ at various Lt,ol 
(left) and T(,oi (right). The binsize is 1/3 dex in both quantities. 
For the observations, only the 112 embedded sources (plotted as 
filled circles on the BLT diagrams) are included. 

shown in Figure [3]). Column 1 of Table [T] lists the Li,oi 
and Ttoi limits of each bin, and column 2 lists the fraction 
of the 112 embedded sources in each bin. For compari- 
son, column 3 lists the fraction of total time the YE05 
model spends in each bin, calculated from Equation 1101 
Columns 4 — 8 list the same thing for the revised models 
presented below. The luminosity problem in the YE05 
model is emphasized by the fact that 76.8% of the ob- 
served sources have 0.1 < Li,oi < 10 Lq while 16.1% have 
10 < Lboi < 1000 L0, whereas the YE05 model spends 
only 17.4% of the time at 0.1 < Lboi < 10 L© but 77.2% 
of the time at 10 < Ltoi < 1000 Lq. 
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Fig. 3. — BLT diagram showing the bins used in the comparison 
between the fraction of total sources and fraction of total time spent 
in each bin presented in Table [T] The thick dashed line shows the 
ZAMS (D'Antona & Mazzitelli 1994) from 0.1 to 2.0 Mq. The 
symbols show the Young Stellar Objects from Evans et al. (2009) 
in this diagram and hold the same meaning as in Figure [Tl 



4. MODIFICATIONS TO THE ORIGINAL MODEL 

We now make several modifications to the YE05 model 
in a step-by-step fashion, adding them one at a time and 
discussing the results of each before adding in the next. 
To be specific, (1) we include isotropic scattering off dust 
grains, (2) we include a circumstellar disk directly in the 
radiative transfer, (3) we include a rotationally flattened 
protostellar envelope density structure, (4) we include 
the effects of mass-loss and outflow cavities, and (5) we 
include episodic accretion. These additions are summa- 



rized in Table [H All of these additions are made possi- 
ble by switching to the two-dimensional RADMC rather 
than the one-dimensional dust radiative transfer package 
DUSTY (Ivczic et al. 1999) used by YE05. 

4.1. Model 1: Scattering 

YE05 assumed the dust opacities of Ossenkopf & Hen- 
ning (1994) appropriate for thin ice mantles after 10^ 
yr of coagulation at a gas density of 10^ cm^'^ (0H5 
dust), which several recent studies have shown to be ap- 
propriate for cold, dense cores (e.g., Evans et al. 2001; 
Shirley et al. 2002; Young et al. 2003; Shirley et al. 2005). 
These opacities do not extend below 1.25 /xm, and they 
include only the total dust opacity (^tot) rather than 
the contributions from absorption (^ahs) and scattering 
[liscat) separately. To remedy this, YE05 also used the 
dust opacities calculated by Pollack et al. (1994) for dust 
grains with a radius of 0.1 /im at a temperature of 10 K, 
which give Kabs and Kgcat separately and extend down to 
0.091 /im. Noting that Ktot according to 0H5 and Pollack 
et al. agreed at short wavelengths, YE05 simply obtained 
Kabs and Kscat from Pollack et al. shortward of 1.25 /im. 
Longward of 1.25 /im, they used Ktot from the 0H5 dust 
and the albedo (ratio of Kscat/ i^tot) from Pollack et al. to 
apportion the total 0H5 opacity among scattering and 
absorption components. 

Unfortunately, YE05 were not able to include the ef- 
fects of scattering when using the 1-D radiative transfer 
code DUSTY to calculate the envelope dust temperature 
profile and final SED of the protostar-|-disk-|-envelope 
system. DUSTY assumes that scattering from dust 
grains is isotropic, when in reality the grains preferen- 
tially forward-scatter light. As described in more detail 
by YE05, assuming isotropic scattering results in an ar- 
tifical near-infrared peak in the SED of any core exposed 
to the ISRF, even a starless core, from backscattering 
of ISRF photons. This peak can be as strong as the 
true peak from thermal dust emission in the far-infrared 
and submillimeter (YE05). As the ISRF is the dominant 
source of heating at early times, YE05 were forced to 
ignore the effects of scattering in order to produce rea- 
sonable results. 

However, by ignoring scattering, they removed an im- 
portant opacity source from their model. Figure|31 which 
plots Kabs, Kscat, and the ratio of Kscat to Kabs as a func- 
tion of wavelength A, shows that the opacity from scat- 
tering dominates that from absorption over the approxi- 
mate wavelength range 0.1 < A < 10 fim. As this is the 
same approximate wavelength range where the emission 
from the protostar-|-disk input SEDs peak, neglecting the 
opacity from scattering will artificially boost the amount 
of short-wavelength radiation escaping from the models. 
Other, more recent studies of individual sources using 
DUSTY to model the observed SEDs have attempted to 
correct for this by treating the total opacity {Kabs + Kscat) 
entirely as absorption (e.g., Bourke et al. 2006; Dunham 
et al. 2006). This method will give the correct amount 
of total opacity, although it will overcorrect and produce 
too little short-wavelength radiation since some of the 
absorbed radiation should have been scattered instead. 

Here we revisit the issue of including scattering in 
the radiative transfer. Even with preferential forward- 
scattering off dust grains, at least some near-infrared 
emission is still expected. Indeed, Foster & Goodman 
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Fig. 4. — Top: Absorption and scattering opacities (/tabs a^nd 
f^scat, respectively) constructed by YE05 based on the dust opacity 
models of Ossenkopf & Henning (1994) and Pollack et al. (1994). 
Bottom: Ratio of scattering to absorption opacity. 

(2006) detect extended near-infrared emission arising 
from such scattering, which they call "cloudshine" , in 
very deep near-infrared images of dark clouds in Perseus. 
Unlike near-infrared emission arising from a protostar, 
which is compact in nature, the cloudshine originates 
from the full extent of the core. Thus, assuming typical- 
sized apertures of a few arcseconds or less are used, only 
a small amount (approximately equal to the ratio of the 
solid angle subtended by the aperture to that subtended 
by the full extent of the core) of this emission would 
be included in near-infrared photometry of embedded 
sources. Furthermore, subtraction of the sky background 
performed in any standard photometry procedure would 
remove the small amount of cloudshine that is included 
in the aperture. Thus, none of the cloudshine should be 
included in these models if they are to be compared to 
observations. 

By switching to the 2-D dust radiative transfer package 
RADMC, we are able to simulate observations by includ- 
ing both small apertures and background subtraction. 
Like DUSTY, RADMC assumes that the scattering pro- 
cess is isotropic. However, as RADMC is a Monte Carlo 
code that follows individual photons from their creation 
at the central source through to their final escape from 
the system, the locations of the source of the observed 
photons are preserved. We thus calculate the final SED 
in apertures of fixed sizes to include the emission from 
the compact source but exclude most of the diffuse emis- 
sion from scattering of the ISRF. Following Crapsi et 
al. (2008), we select the aperture radii to approximately 
match the resolution of the Spitzer Space Telescope: 2" 
(280 AU at 140 pc) for A < 10 /im, 6" (840 AU at 140 
pc) for 10 < A < 40 Aim, and 20" (2800 AU at 140 pc) 
for 40 < A < 100 /im. Longward of 100 ^m we use an 
aperture large enough to encompass the full extent of the 
model. To simulate removing any remaining cloudshine 
through sky background subtraction, we run the entire 
model grid a second time, including only the ISRF and 
no internal source of luminosity. To construct the final 
SED at each timestep, we then subtract the model with 
only external heating from the model with both inter- 




Wavelength (microns) 

Fig. 5. — Model 1 spectral energy distributions (SEDs) at select 
ages through the 224,000 yr collapse o f the 1 Mq core. In each 
panel, the solid line shows the model 1 ( i|4.1l l SED while the dotted 
line shows the original YE05 SED. There is no YE05 SED shown 
in the last panel because the YE05 1 Mq model finished collapsing 
at 210,000 yr. 

nal and external heating for A < 10 /xm, since no core 
heated only externally by the ISRF will emit any signif- 
icant thermal radiation at such short wavelengths (e.g., 
Evans et al. 2001). 

Figure [5] compares the SEDs with scattering included 
to the YE05 SEDs at various times throughout the col- 
lapse of the 1 M0 core (analogous to Figure 8 in YE05). 
As expected, the SEDs with scattering included have sig- 
nificantly less short-wavelength emission than found by 
YE05. At later times there is also an increasing discrep- 
ancy in the long- wavelength emission. Since the emission 
at such wavelengths is optically thin it traces the total 
envelope mass, signifying a growing difference with time 
in the total remaining envelope mass. This is caused by 
the 5% lower Menv in this work compared to the YE05 
model (see pTTjl . 

Figure [6] shows the evolution of Lboi, Tboi, and 
Lboi/ Lsmm, both for the original YE05 model and for 
model 1, plotted against the ratio of internal (pro- 
tostar -|-disk) to total (protostar -|-disk-|-envelope) mass 
(Mint/Mtot)- Based on the physical definition of Stage 
as the portion of the embedded phase when at least 
half of the total system mass is still in the envelope, the 
model should cross the Class O/I boundaries in Tboi and 
Lboi/Lsmm when M^nt/Mtot =0.5. From this, YE05 
concluded that Lboi/Lsmm is a much more reliable evo- 
lutionary indicator than Tboi- Indeed, Figure [6] shows 
that, in the YE05 model, the T^oi Class O/I boundary is 
crossed when Mmt/Mtot = 0.48, 0.15, and 0.05 for the 
0.3, 1, and 3 Mq cores, respectively. On the other hand, 
the YE05 model crosses the Lboi/Lgmm Class O/I bound- 
ary when Mint/Mtot = 0.60, 0.51, and 0.40, respectively. 

Including the opacity from scattering changes these re- 
sults. As noted above, including scattering significantly 
decreases the short-wavelength emission since the opac- 
ity at these wavelengths is increased by up to a factor 
of 10 (Figured]). As a consequence, the calculated Tboi 
of a given model decreases substantially. To be quan- 
titative, for the 1 Mq core, except for very early times 
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Fig. 6. — Observational signatures Ljoj (top panels), T^oj 
(middle panels), and I'holl^smm (bottom panels) versus 
Mint/Mtot, the ratio of internal (protostar+disk) to total (proto- 
star+disk+envelope) mass. The left panels show the YE05 model 
results while the right panels show the model 1 results. The class 
boundaries in Tjoj are taken from Chen et al. (1995) while the class 
divisions in Li,oi/Lamm are taken from YE05. The discontinuities 
in both Tjoj and Li,ai/Lsmm are artifacts introduced by the switch 
from the Shu (1977) density profile to a power-law density profile. 
They are present with the same magnitude in the YE05 model but 
are less readily apparent because of the logarithmic axis scale. 



when neither model features short-wavelength emission, 
the model with scattering included has a calculated Ti,oi 
lower by a factor of about 1.5 — 6. This reduction in Ti,oi 
is evident in Figure [SI For the model with scattering 
included, the Tboi Class O/I boundary is crossed when 
M,nt/Mtot = 0.70, 0.28, and 0.09 for the 0.3, 1, and 3 
M© cores, respectively, uniformly later than in the YE05 
models. This same model crosses the Lboi/ Lsmm Class 
O/I boundary when M.^t/Mtot = 0.66, 0.54, and 0.22, 
respectively. 

Based on these results, we conclude that it is im- 
portant to include the opacity from scattering at short 
wavelengths. Doing so reduces the amount of short- 
wavelength emission and is thus crucial for comparing 
model and observed SEDs. It also reduces the calculated 
Tboi of a given model by a factor of '--^ 1.5 — 6 depend- 
ing on the exact parameters of the model. It does not 
significantly affect iboz/^smm, however, since this quan- 
tity is much less sensitive to the exact amount of short- 
wavelength emission (the lower values of Lboi/Lsmm at 
late times compared to the YE05 model is a consequence 
of the increased long-wavelength emission arising because 
of the 5% lower mass accretion rate, as described above, 
and is unrelated to including the opacity from scatter- 
ing). The Class O/I Tfjoj boundary is still crossed too 
early for the 1 and 3 M© cores compared to the Stage 
/I boundary, but the discrepancy is not as bad as in the 
YE05 model. While Lboi/ Lsmm remains the most reli- 
able evolutionary indicator for associating physical Stage 
with observational Class, we caution that, in practice, it 
is more difficult to calculate and significantly more prone 
to error than Thoi, depending on the exact submillimeter 
wavelengths at which observations are available (Dun- 
ham et al. 2008). 

Figure [7] shows a BLT diagram for model 1, similar to 
Figured] for the YE05 model. While including the opac- 
ity from scattering has important consequences, as de- 
scribed above, the general luminosity problem described 
in fJT] remains. 

Figure [5] shows Lf,oi and Thoi histograms for both the 
observations and model 1 , similar to Figure [2] for the 
YE05 model, and column 4 of Table [T] lists the fraction 
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Fig. 7. — Same as Figure [T] except now for model 1 rather than 
the YE05 model, and without model tracks. The grayscale pixels 
indicate the fraction of total time the model spends in each L^o; — 
Tboj bin, calculated from Eouation llOl The grayscale is displayed in 
a logarithmic stretch with the scaling chosen to emphasize the full 
extent of the models in Li,ai — Ti,ai space. The mapping between 
grayscale and fraction of total time is indicated in the legend. The 
class boundaries in Tfj^i are taken from Chen et al. (1995). The 
thick dashed line shows the ZAMS (D'Antona & Mazzitelli 1994) 
from 0.1 to 2.0 M©. The colored symbols show the Young Stellar 
Objects from Evans et al. (2009) in this diagram; the colors and 
symbols hold the same meaning as in Figure [l] 

of total time model 1 spends in various Lboi — Tboi bins. 
These results again emphasizes the luminosity problem. 
76.8% of the observed sources have 0.1 < Lboi < 10 L© 
while 16.1% have 10 < Lboi < 1000 L©, whereas the 
models spends only 13.4% of the time at 0.1 < Lboi < 10 
L© but 80.6% of the time at 10 < Lboi < 1000 L©. Fur- 
thermore, a K-S test shows that there is less than a 0.1% 
probability that the observed and model Lboi histograms 
represent the same underlying distribution. A similar 
K-S test gives a 42% probability that the observed and 
model Tboi histograms represent the same underlying dis- 
tribution. The increase in short-wavelength opacity and 
corresponding decrease in both short-wavelength model 
emission and model Tboi is clearly seen in the Tboi his- 
togram in that model 1 spends more time at low Tboi 
{Tboi ^5 200 K) compared to the YE05 model. Compared 
to the observations, model 1 overpredicts the fraction of 
sources observed at these low values of Tboi ■ 
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Fig. 8. — Same as Figure [2] except now for model 1 rather than 
the YE05 model: histograms showing the fraction of total sources 
(observations; solid lines) and fraction of total time spent (models; 
dashed lines; calculated from Equation I10| l at various Ltot (left) 
and Tjoj (right). The binsize is 1/3 dex in both quantities. For 
the observations, only the 112 embedded sources (plotted as filled 
circles on the BLT diagrams) are included. 



4.2. Model 2: Two-Dimensional Disk 

A protostellar disk is inherently a two (or higher) di- 
mensional object, but the radiative transfer in the YE05 
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model was calculated with DUSTY, a one-dimensional 
radiative transfer code. In their model and model 1 
above, a disk was included by calculating the emission 
from a disk with specified parameters, averaging this 
emission over all inclinations, and then adding the re- 
sult to the (proto)stellar spectrum to form the final in- 
put SED for the radiative transfer calculation (Butner et 
al. 1994; Adams et al. 1988). However, since this work is 
performed with RADMC, a 2-D radiative transfer code, 
we are able to include a disk directly in the radiative 
transfer rather than simulate its presence by the above 
method. 

Several other recent authors have published 2-D mod- 
els of embedded protostars that include disks (e.g., Whit- 
ney et al. 2003a; Whitney et al. 2003b; Robitaille et 
al. 2006; Crapsi et al. 2008). Following their examples, 
we include a relatively simple disk with a density pro- 
file featuring a power-law in the radial coordinate and a 
Gaussian in the vertical coordinate: 



Pdiskis, z) = Pa 



exp 



- - 

2 \H, 



(11) 



In equation 1111 z is the distance above the midplane 
(z = rcos6, with r and 9 the usual radial and zenith angle 
spherical coordinates) and s is the distance in the mid- 
plane from the origin (s — ^Jr^ — z'^). The quantity 

is the disk scale height and is given by Hs = Ho (t") ; 

where i?o is the scale height at the reference midplane 
distance So- Following Whitney et al. (2003a), we set 
Eo = 10 AU at So = 100 AU. The parameter /3 de- 
scribes how the scale height changes with s and sets the 
flaring of the disk, while a describes how the midplane 
density profile varies with s. Again following Whitney 
et al. (2003a), we choose /3 = 1.25 and a = 2.25. These 
values are close to those adopted by Crapsi et al. (2008; 
/3 = 9/7, a = 16/7) to correspond to the self-irradiated 
passive disk model of Chiang & Goldreich (1997), and 
also to the best-fit values found by Sauter et al. (2009; 
/3 = 1.4, a — 2.2) from a detailed model of the edge-on 
circumstellar disk CB 26. The mass of the disk is con- 
trolled by the parameter po, where the mass, inner disk 
radius, and outer disk radius evolve following the descrip- 
tion in ii2.1l The disk surface density profile, calculated 
by integrating Equation [Tl] over the vertical coordinate 
has a radial power-law index of S(s) oc s~^, where 
p = a — (3. With our adopted values of a and /3, p = 1. 

We include scattering as described in ii4.ll and we in- 
clude the same luminosity components described in ^2.'2^ 
In the original model and that presented in ! j4.11 LeioD, 
LdioPi and L^m are treated as intrinsic disk luminosity, 
while LeioP, ^_Di?,, and Lphot are treated as protostellar 
luminosity. However, our 2-D radiative transfer pack- 
age RADMC is limited to one internal source of pho- 
tons (the protostar) and one external source (the ISRF). 
The disk in this model is thus treated as a purely repro- 
cessing disk; it reprocesses radiation from the protostar 
but does not have its own intrinsic luminosity. To keep 
the overall model luminosities correct, we treat all six 
sources of luminosity as protostellar and calculate the 
resulting temperature profiles and emergent SEDs of the 
disk+envelope model. The main consequence of adopt- 
ing this purely reprocessing disk is less mid-infrared emis- 




Wavelength (microns) 

Fig. 9. — Same as Figure |51 except now for model 2: SEDs at 
select ages through the collapse of the 1 Mq core. In each panel, 
the lines show the model 2 SEDs at various inclinations, with the 
line weight key given in the first panel. The dotted line shows, for 
comparison, the SED at that time for model 1. 

sion, as discussed in greater detail below. 

Before presenting our results, we note that the equa- 
tions presented by Adams & Shu (1986) for the luminos- 
ity components assume a spatially thin disk (confined 
to the z — plane), which is clearly not the case for 
the fiared disk adopted here. Material accreting onto 
the disk will not fall all the way to the midplane be- 
fore joining the disk and thus will not fall as deep into 
the potential well. As a consequence, LEtoD, calculated 
assuming a thin disk, is actually an upper limit to the 
true value. However, the total model luminosity should 
be correct, since the decreased luminosity from accretion 
onto the disk should be compensated by the increased 
luminosity from the accretion of this material from the 
disk onto the protostar. The YE05 model also featured 
this small physical inconsistency since the parameters of 
their simulated, 1-D disk were chosen to simulate a fiared 
disk. The overall effects on the results should be negli- 
gible since we are interested in the global evolution of a 
collapsing protostellar core rather than the details of the 
individual components of the total luminosity. 

Figure [9] compares the model SEDs at various inclina- 
tions for the 1 M0 core to the model 1 SEDs. There 
is very little inclination dependence at early times when 
the disk is small and not very massive. At late times, as 
the disk grows and becomes more massive, the inclina- 
tion dependence increases, although the only substantial 
change with inclination occurs for nearly edge-on lines of 
sight that pass through the disk. 

Compared to the SEDs with a 1-D, simulated disk, 
the main difference of this model (aside from the inclina- 
tion dependence) is a decrease in the 3—10 /im emis- 
sion. This difference arises because model 1 includes 
intrinsic disk luminosity which raises the disk temper- 
ature profile above that arising solely from reprocessing 
of (proto)stellar radiation, whereas model 2 presented 
here is limited to a reprocessing disk only. While the to- 
tal internal luminosity of the model remains correct, as 
described above, there is less mid-infrared emission be- 
cause the disk is generally cooler when heated only by 
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Fig. 10. — Same as Figure[6l except now with model 1 in the left 
panels and model 2 in the right panels: observational signatures 
^bol (top panels), T^oj (middle panels), and Lf,al/Lsmm (bottom 
panels) versus Mint/Mtot, the ratio of internal (protostar+disk) to 
total (protostar+disk+envelope) mass. The different colors show 
the results for different inclinations; the color key is given in the 
upper right panel. The class boundaries in Tjoj are taken from 
Chen et al. (1995) while the class divisions in Lij^i/ Lsmm are taken 
from YE05. The discontinuities in both T(,oj and Li,oi/L 

smm are 

artifacts introduced by the switch from the Shu (1977) density 
profile to a power-law density profile. 

reprocessing. However, as shown below, this difTerence 
has only a small effect on calculated evolutionary indi- 
cators and is not important in the context of this work. 
The exact amount of mid-infrared emission in the final 
SEDs can be adjusted by varying the degree of flaring of 
the disk (through the parameter (3), which will affect the 
amount of (proto)starlight intercepted and reprocessed 
by, and thus the temperature profile of, the disk. 

Figure [To] shows the observational signatures L&oi, Tfco;, 
and Lboi/Lsmm plotted against the ratio of Mmt/Mtot, 
with different colors corresponding to different inclina- 
tions. As expected from the SEDs, there is essentially no 
inclination dependence visible in any quantity except at 
late times, and even then most of the dependence is seen 
in the nearly edge-on line-of-sight (i = 75°) relative to 
the other lines of sight. To be quantitative, for the 1 
core, Thoi calculated from SEDs at the same time viewed 
at 5° and 85° vary by < 5% for all times up to -150,000 
yr, and by < 25% for the remaining times. The ratio of 
Lboi/Lsmm shows similar results: Lboi/Lsmm calculated 
from SEDs viewed at 5° and 85° varies by < 25% except 
at very late times (approximately the last 10^ yr), where 
the nearly edge-on lines-of-sight passes through a disk 
that has become so optically thick that the calculated 
Lboi significantly decreases (evident in the last panel of 
Figure [H]). At these late times, Lboi/Lsmm changes by 
about a factor of four from pole-on to edge-on inclina- 
tions. 

Compared to model 1, the addition of a 2-D disk in 
the radiative transfer causes small reductions in both 
Tboi and Lboi/Lsmm- This is caused by the change to 
a reprocessing disk described above, which shifts some 
of the shorter-wavelength, 3—10 /im emission to longer 
wavelengths and thus decreases both evolutionary indica- 
tors. The effect is small; at any given time, Tboi decreases 
by < 25% and Lboi/Lsmm decreases by < 15%. As a 
consequence, model 2 crosses the Class O/I boundary^*^ 

Different inclinations cross the Class O/I boundary at slightly 
different times, and thus at slightly different values of Mint/Mtot- 
The values presented here are calculated by taking a weighted av- 
erage of the values at each inclination, with the inclination weights 



slightly later: the Tboi Class O/I boundary is crossed 
when Mint/Mtot = 0.76, 0.31, and 0.11 for the 0.3, 1, 
and 3 cores, respectively, and the Lboi/Lsmm Class 
O/I boundary is crossed when Mint/Mtot = 0.74, 0.61, 
and 0.23, respectively. These times are slightly later 
than model 1 (see M.ip , but the overall conclusions about 
the connection between physical Stage and observational 
Class, as defined by Tboi and Lboi/Lsmm, remain un- 
changed. 

A careful inspection of Figure \W\ reveals small-scale 
oscillations in both Tboi and Lboi/Lsmm at late times. 
These are artifacts introduced by the model gridding. 
The grid is logarithmically spaced in radius to ensure 
that there are enough gridpoints at small radii where the 
optical depth is large. As a consequence, the spacing 
between gridpoints becomes relatively large ('^ 5 — 10 
AU per gridpoint) at 50 — 100 AU, the range of disk 
outer radii at these late times. As the disk radius grows 
with time, it generally takes a few timesteps before it 
increases enough to "jump" to the next gridpoint. In the 
model it thus remains fixed at the previous gridpoint, 
but as the mass is increasing at each timestep, the opti- 
cal depth through the disk is also increasing. Once the 
radius increases enough to jump to the next gridpoint, 
the optical depth suddenly decreases slightly, affecting 
the temperature profile, emergent SEDs, and thus cal- 
culated evolutionary indicators. These effects are small 
enough to have a negligible impact on the results. 

Figures [Tl] and [T2] show a BLT diagram and Lboi and 
Tboi histograms for model 2, respectively, and column 5 
of Table [T] lists the fraction of total time model 2 spends 
in various Lboi — Tboi bins. The inclination dependence 
introduced by the disk results in slightly more spread in 
model coverage in Lboi — Tboi space, and shifts the peak 
of the model luminosity distribution to slightly lower lu- 
minosities, but the main model 1 conclusions that the 
model overpredicts both the time spent at high luminosi- 
ties (> 1 Lq) and the time spent at low Tboi (^ 200 K) 
remain unchanged. Indeed, K-S tests on Figure [12] give 
similar results as those performed on Figure [8] there is 
again less than a 0.1% probability that the observed and 
model Lboi histograms represent the same underlying dis- 
tribution, and a 34% probability of the same for the Tboi 
histograms. 

4.3. Model 3: Two- Dimensional Envelope 

All models presented up to this point have featured the 
spherically symmetric density profiles calculated by Shu 
(1977) for the collapse of singular isothermal spheres. 
In reality, a collapsing core that is initially spherically 
symmetric will not remain so if the core has any initial 
rotation; conservation of angular momentum will create a 
rotationally flattened structure. As the models presented 
here do feature initial rotation (set by the initial angular 
velocity fip; see ii2.1[) . this effect should be included. The 
move to the 2-D radiative transfer code RADMC enables 
us to do this. 

To include the effects of rotation on the evolution of 
the model, we adopt the solution for the collapse of a 
slowly rotating core presented by TSC84. The core is ini- 
tially a spherically symmetric, singular isothermal sphere 
with a density distribution n cx r~^^ , identical to the Shu 

as defined in the discussion following Equation llOl ( ^4.111 . 
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Fig. 11. — Same as Figure [T] except now for model 2. The 
grayscale pixels indicate the fraction of total time the model spends 
in each L),oi — Ti,oi bin, calculated from Eguation llOl The grayscale 
is displayed in a logarithmic stretch with the scaling chosen to em- 
phasize the full extent of the models in Liy^i — Tf,ai space. The 
mapping between grayscale and fraction of total time is indicated 
in the legend. The class boundaries in Tjoj are taken from Chen 
et al. (1995). The thick dashed Une shows the ZAMS (D'Antona 
& MazziteUi 1994) from 0.1 to 2.0 Mq. The colored symbols show 
the Young Stellar Objects from Evans et al. (2009) in this diagram; 
the colors and symbols hold the same meaning as in Figure [7] 
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Fig. 12. — Same as FigurefS] except now for model 2: histograms 
showing the fraction of total sources (observations; solid lines) and 
fraction o f to tal time spent (models; dashed lines; calculated from 
Equation IIOI I at various Li,o! (loft) and Tj,oj (right). The binsize 
is 1/3 dex in both quantities. For the observations, only the 112 
embedded sources (plotted as filled circles on the BLT diagrams) 
are included. 

(1977) solution. The outer radius is again truncated to 
set the initial core mass. As collapse proceeds, the solu- 
tion takes on two forms: an outer solution that is similar 
to the non-rotating, spherically symmetric solution and 
an inner solution that exhibits flattening of the density 
profile. Since material falling in to the central regions 
originates from larger radii and thus carries more an- 
gular momentum as time progresses, the radius where 
the inner solution must be used, and thus the radius at 
which flattening becomes significant, increases with time 
{rfiat oc f^gt^; TSC84), reaching - 2000 AU at late times 
for all three initial mass cores (the longer collapse times 
of the higher mass cores are offset by slower initial angu- 
lar velocities). 

Once half of the initial core mass has accreted onto the 
protostar-|-disk system, the infall radius, which moves 
outward at the sound speed, exceeds the envelope outer 
radius. If we simply continued to use the TSC84 solution, 
"extra mass" that originated beyond the outer radius 
would begin to collapse and eventually move within this 
radius. Thus, there would still be mass remaining in the 
envelope once the initial mass of the core has accreted, 
which is clearly not self-consistent (although see Myers 
[2008] for an analytic model that includes protostellar ac- 



cretion from a core embedded in a uniform background 
that also partially accretes onto the protostar). YE05 
also faced this problem in their 1-D model and avoided 

it by simply switching to an n oc Tenv power-law density 
profile once the infall radius exceeded the envelope outer 
radius, with the outer radius kept fixed and the desired 
mass at each timestep in the model evolution achieved 
by adjusting the normalization of the power-law. How- 
ever, this results in a model that is no longer physically 
self-consistent, since it does not make sense to have a 
core that is collapsing at all radii (which is the case once 
the infall radius exceeds the outer radius) and thus de- 
creasing in mass but remaining a fixed size. Furthermore, 
such a solution is not available to us here since we wish to 
retain the feature of the TSC84 solution that the radius 
at which the density profile exhibits fiattening increases 
with time. Thus, as an alternative, we use the veloc- 
ity profiles given by the TSC84 solution and allow the 
outer radius to decrease once the infall radius exceeds 
the initial outer radius. This is illustrated by Figure [T3l 
which shows the outer radius of the core and inward ra- 
dial velocity at this outer radius as a function of time 
for the three different initial mass cores. The effect of 
this change in the calculation of the density profiles for 
the second half of the collapse of each core is discussed 
below. 
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Fig. 13. — The model 3 envelope outer radius (top panel) and 
inward radial velocity at this outer radius (bottom panel) as a 
function of time for the 0.3, 1, and 3 Mq cores. The outer radius 
of each core remains constant with no inward radial velocity for 
the first half of the collapse when the infall radius has not yet 
reached the outer boundary of the core, and decreases with an 
increasing inward radial velocity once the infall radius reaches the 
outer boundary. 



Figure [HI compares the model 3 SEDs at various incli- 
nations for the 1 Mq core to the model 1 SEDs presented 
in i 34.ll As was the case for model 2, the inclination de- 
pendence is small at early times and increases throughout 
the collapse of the core. SEDs at late times show notice- 
ably less short-wavelength emission. Part of this is due 
to the change to a purely reprocessing disk, as discussed 
above. However, this deficit in short-wavelength emission 
becomes especially noticeable at very late times and has 
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a second cause: the decrease in the outer radius of the 
envelope as material at the initial outer edge collapses. 
The reason for this is best understood by assuming this 
is a 1-D, spherically symmetric model with a power-law 
density profile given by p{r) = pf{r/rf)~P, with pf the 
density at a fiducial radius rf. Assuming that 1 < p < 3, 
the equations for envelope mass and optical depth are: 




Wavelength (microns) 

Fig. 14. — Same as Figure [9] except now for model 3: SEDs 
at select ages through the collapse of the 1 Mq core. In each 
panel, the lines show the model 3 (scattering, disk, and rotationally 
flattened envelope included) SEDs at various inclinations, with the 
line weight key given in the first panel. The dotted line shows, for 
comparison, the SED at that time for model 1. 
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In comparison to model 1 at the same timestep, the 
envelope has a smaller outer radius but identical mass, 
which, according to Equation 1121 requires a larger value 
of Pf, the normalization of the power-law. According 
to Equation I13( since the envelope inner radius is un- 
changed, this increases the optical depth through the en- 
velope, giving less short-wavelength emission. The situ- 
ation is slightly more complicated in reality since we are 
comparing a 1-D, spherically symmetric model to a 2-D, 
rotationally flattened model, but the analogy holds. Al- 
lowing the outer radius to decrease following the TSC84 
collapse solution increases the optical depth compared 
to YE05 and models 1 — 2 at the same timestep, where 
the outer radius is held fixed, and causes more of the 
short-wavelength emission to be reprocessed to longer 
wavelengths than for these previous models. 

Figure [15] shows the observational signatures Ljjoi, Thoi, 
and Lboi/Lsrnm plotted against the ratio of Mint/Mtot, 
with different colors corresponding to different inclina- 
tions. Although neither Tboi nor Lboi/ Lsmm increases 
to values as large as those reached in model 2 at late 
times due to the decrease in short-wavelength emission 
described above, most of the evolution is similar to 



model 2. The Tboi Class O/I boundary is crossed^ ^ when 
A'hnt/Mtot = 0.68, 0.33, and 0.12 for the 0.3, 1, and 
3 M0 cores, respectively, and the Lhoi/Lsmm Class O/I 
boundary is crossed when Mtnt/Mtot = 0.73, 0.61, and 
0.24, respectively. Comparison to the model 2 results 
shows that the overall conclusions about the connection 
between physical Stage and observational Class, as de- 
fined by Tboi and Lboi/Lsmm, again remain unchanged. 
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Fig. 15. — Same as Figure [TOl except now for model 3: obser- 
vational signatures L^oj (top panels), T^qj (middle panels), and 
Lbol/Lsmm (bottom panels) versus Mint / Mtot, the ratio of in- 
ternal (protostar+disk) to total (protostar-|-disk-|-envelope) mass. 
The left panels show the model 1 results while the right panels show 
the model 3 results. The different colors show the results for dif- 
ferent inclinations; the color key is given in the upper right panel. 
The class boundaries in T^oj are taken from Chen et al. (1995) 
while the class divisions in Lf,ai/Lsmm are taken from YE05. The 
discontinuities in both T^oi and Li,oi / Lsmm for model 1 are arti- 
facts introduced by the switch from the Shu (1977) density profile 
to a power-law density profile. 



Both quantities show an inclination dependence that 
increases with time. Using the 1 Mq core as an example, 
Tboi calculated from SEDs viewed at 5° and 85° vary by 
< 5% only for the first 96,000 yr (compared to the first 
150,000 yr for model 2), and the variation increases to 
55% for very late times in the model evolution (compared 
to 25% for model 2). This inclination dependence is not 
confined solely to lines-of-sight viewed close to edge-on 
since the flattening of the envelope affects all viewing 
angles: Tboi calculated from SEDs viewed at 5° and 25° 
vary by up to 20% at late times, and Tboi calculated 
from SEDs viewed at 5° and 45° vary by up to 40% 
at late times. For comparison, model 2 showed < 5% 
variations in Tboi at any given time for all inclinations < 
70° . Similar results are found for the ratio of bolometric 
to submillimeter luminosity. Lboi / L g„^m calculated from 
SEDs viewed at 5° and 85° vary by up to 50% at late 
times, compared to 25% (except for the last 10,000 years) 
for model 2. 

While the inclusion of a disk and rotationally flattened 
envelope following TSC84 does induce a moderate incli- 
nation dependence, we note here that it is a much smaller 
dependence than found by other authors investigating 
2-D models of embedded protostars (e.g., Whitney et 
al. 2003a; Whitney et al. 2003b; Robitaille et al. 2006; 
Crapsi et al. 2008). For example, Crapsi et al. (2008) 
found a variation in Tboi with inclination ranging from 
factors of ~ 2 — 5 depending on the exact model param- 
eters. Most of the explanation for this difference resides 

Again calculated as a w eighted average of the different incli- 
nations, as described in i|4.2l 
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in the fact that we do not yet include outflow cavities, 
whereas other authors do. This will be addressed in m.4l 
However, there is a second, important difference in our 
models that is worth pointing out: we follow the exact 
TSC84 solution whereas these other studies do not. In- 
stead, they adopt density profiles that exhibit rotational 
flattening at all radii. As noted by Terebey et al. (2006), 
these other models only agree with the TSC84 solution 
at small radii where the inner solution is valid; at large 
radii they diverge. 

Figures [H] and [T7] show a BLT diagram and Lboi and 
Tboi histograms for model 3, and column 6 of Table [1] lists 
the fraction of total time model 3 spends in various Lboi 
— Tboi bins. The decrease in T^oi at late times due to al- 
lowing the outer radius of the core to shrink as material 
initially at this outer radius collapses is seen in Figure [161 
in that the models do not extend beyond about 1000 K, 
and in Figure [17] in the increased amount of time spent 
at Tboi < 200 K relative to that spent at Tboi > 200 K. 
However, the luminosity distribution remains essentially 
unchanged, and the main model 1 and model 2 conclu- 
sions that the model overpredicts both the time spent at 
high luminosities (> 1 Lq) and the time spent at low 
Tboi 200 K) remain unchanged. Indeed, K-S tests 
on Figure [T7| give similar results as those performed on 
Figure [T^l there is again less than a 0.1% probability 
that the observed and model Lboi histograms represent 
the same underlying distribution, and a 13% probability 
for the Tboi histograms. The lower probability that the 
Tboi distributions are the same compared to models 1 — 2 
(13% versus ^ 35 — 40%) arises from the shift to even 
more time spent at low Tboi in the model compared to 
observations. 
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Fig. 16. — Same as Figure 1111 except now for model 3. The 
grayscale pixels indicate the fraction of total time the model spends 
in each L^oj — T;,;,; bin, calculated from Equation llOl The grayscale 
is displayed in a logarithmic stretch with the scaling chosen to 
emphasize the full extent of the models in L;,;,; — r^oj space. The 
mapping between grayscale and fraction of total time is indicated 
in the legend. The class boundaries in Tjjo; are taken from Chen 
et al. (1995). The thick dashed Une shows the ZAMS (D'Antona 
& MazziteUi 1994) from 0.1 to 2.0 Mq. The colored symbols show 
the Young Stellar Objects from Evans et al. (2009) in this diagram; 
the colors and symbols hold the same meaning as in Figure [7] 



4.4. Model 4- Mass-Loss and Outflow Cavities 

All of the models considered so far have featured a 
100% collapse efficiency^^ ; all material initially in the 

We refrain from using the term star formation efficiency here, 
as this term is more commonly used to refer to the fraction of mass 
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Fig. 17. — Same as Figure [T2l except now for model 3: histograms 
showing the fraction of total sources (observations; solid lines) and 
fraction o f to tal time spent (models; dashed lines; calculated from 
Equation I10| l at various L),oi (left) and T(,oi (right). The binsize 
is 1/3 dex in both quantities. For the observations, only the 112 
embedded sources (plotted as filled circles on the BLT diagrams) 
are included. 

core is in the star-|-disk system at the end of collapse. In 
reality, however, this is not the case. Some of the mass is 
instead ejected from the system, entraining and removing 
envelope material as it propagates through the core (e.g., 
Bontemps et al. 1996; Bally et al. 2007; Arce et al. 2007). 
Here we incorporate a simple, idealized process of mass- 
loss and the opening of outflow cavities into the evolu- 
tionary model to test whether or not the increased in- 
clination dependence introduced by outflow cavities can 
add enough spread to calculated bolometric luminosities 
and temperatures to bring the model in agreement with 
observations without resorting to episodic accretion. 

A fraction of the material accreted from the disk to 
the protostar is ejected in the form of a jet or wind. 
The exact value of this fraction, Mj/Macc, depends on 
uncertain models of how the jet /wind is launched. As 
compiled by Bontemps et al. (1996), possible values are 
0.1 - 0.5 (Shu et al. 1994), - 0.1 (Pelletier & Pudritz 
1992), and 10"^ - 0.1 (Wardle & Konigl 1993). More re- 
cently, observations of the variation of jet diameters with 
distance from their driving sources have been consistent 
with models giving Mj/Macc > 0.03 (Ray et al. 2007 and 
references therein). Given the above range, we assume 

Mj/Macc = 0.1; 10% of the mass accreting from the disk 
to the protostar is instead ejected from the system in a 
jet /wind, decreasing the growth of the protostellar mass 
accordingly. 

The jet/wind entrains and removes material as it prop- 
agates through the envelope, driving a molecular outflow 
(e.g., Bontemps et al. 1996; Bally et al. 2007; Arce et 
al. 2007; although see Machida et al. 2008 for an al- 
ternate explanation of the origin of molecular outflows). 
Conservation of momentum gives 



fMjvj 



(14) 



where Mj and Mo are the mass- loss rates of the jet/ wind 
and outflow, respectively, vj and Vo are the jet /wind and 
outflow velocities, respectively, and / is the efficiency 
with which the jet /wind transfers its momentum to the 
ambient medium. We assume a typical jet velocity of 
150 km s~i (Bontemps et al. 1996; Bally et al. 2007), 
consistent with being greater than the 6 — 60 km s~^ 

initially in the core that ends up in the star. By that definition, the 
YE05 model and models 1-3 in this work all featu re st ar formation 
efficiencies of 75% (set by the parameter r)£i; see i|2.1l l. 
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escape velocities of 0.1 — 1 M© protostars assuming jet 
launching radii of 0.5 - 5 AU (Ray et al. 2007). We as- 
sume an outflow velocity of 10 km s~^, consistent with 
observations of outflows with typical velocities of 4 — 5 
km s~^ for low-luminosity Class sources (M. Dunham 
et al. 2009, in preparation; J.E. Lee et al. 2009, in prepa- 
ration), 5—15 km s~^ for embedded sources in Perseus 
(HatcheU et al. 2007, 2009), and -10 km s"! for the 
sample of 45 embedded sources compiled by Bontemps et 
al. (1996). The efficiency with which the jet/wind trans- 
fers its momentum to the ambient medium is not well 
characterized and likely varies with environment (e.g., 
Moraghan et al. 2008). We assume an efflciency of / = 1 
to maximize the effects of mass loss and opening of out- 
flow cavities on the results. 

Given the above assumptions, the amount of enve- 
lope mass entrained by the jet/wind is Mentrained = 
15Mejected- Thus, 10% of the mass accreting from the 
disk to the protostar is instead ejected, and 15 times the 
mass of this ejected material is entrained in the outflow. 
The removal of the entrained material is implemented 
into the model by opening outflow cavities (assumed to 
be completely devoid of material) in the envelope. Fol- 
lowing Crapsi et al. (2008), we assume that the outflow 
cavities follow streamlines of the collapse solution, giving 
funnel-shaped cavities that are conical at large radii. The 
size of the outflow cavity, defined by Oc, the semi-opening 
angle of the cone at large radii (see Crapsi et al. [2008] 
for details), increases with time to remove the mass en- 
trained at each timestep. Assuming spherical symmetry 
of the envelope^'^, the ratio of mass entrained and re- 
moved to the total envelope mass is given by the ratio of 
the solid angle of the cavity opened to the total solid an- 
gle of the envelope, which, assuming a pre-existing cavity 
with size 9c,oid is already present, is given as: 

)if , . , O f^^ f!!''''^°"' sin6 dO d(j) 

^^-1- entrained 2 ^ ^cavity 2 ^ '^^c^old 

Mtotai ^totai J^'^ sin 9 dO dcj) 



{cos Ocold - COS Ocnew) 
(1 -I- cos Uc,old) 

The factor of 2 accounts for the bipolar nature of the 
cavities. The semi-opening angle of the outflow cavity, 
dc,new, is calculated from Equation [15] at each timestep, 
with the entrained mass at that timestep calculated as 
described above and Mtotai the total envelope mass at 
that timestep. 

The opening of outflow cavities causes a decrease in 
Mf^nv) the rate at which material is accreting from the 
envelope onto the protostar -|-disk system, since material 
is no longer accreting from the full 47r steradian. As 
above, the mass accretion rate is calculated from the ratio 
of solid angles as follows: 

While the envelope density profiles used here are t he TSC84 
profiles and are thus not spherically symmetric (see i|4.3| l. spherical 
symmetry is still a good approximation at large radii where most 
of the mass is located. 
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Fig. 18. — Evolution of the masses of various model components 
(top panel), outfiow cavity semi-opening angle (middle panel), and 
envelope mass accretion rate (bottom panel) versus time for the 
model 4 1 Mq initial mass co re, w hich includes mass-loss and out- 
fiow cavities as described in i|4.4l 



M'en^ , o ^ca.Uy ^ Jo^C ^^^l 9 dO d4> 

—. — = i — / — — = i — Z — 7^ = COS Ur 

MZlf hiatal J^^ /; sm 9 d9 d<j) 

(16) 

where the factor of 2 again accounts for the bipolar na- 
ture of the cavities and M°^j9 = 4.57x 10"*^ Mq yr^^. 
We neglect the effects of opening cavities and remov- 
ing mass on the collapse solution itself; the overall ef- 
fect would be to slow the collapse of the core (see Myers 
2008). 

Figure [18] shows the effects of including mass loss and 
outflow cavities as described above for the 1 Mq initial 
mass core. Plotted are the masses of the various model 
components (protostar, disk, envelope, and outflow), 9c, 
and Menv versus time. The envelope mass decreases with 
time, with a change in slope to a faster decrease once the 
disk forms and mass is ejected. At this stage the outflow 
mass begins to grow, reaching 0.4 Mq by the end of col- 
lapse. The envelope mass accretion rate decreases as the 
outflow cavity increases in size. Eventually 9c reaches 
90° and the collapse ends at 165,000 yr. Compared to 
the 224,000 yr collapse duration when no mass-loss is 
included, this model is —25% shorter in total duration. 
Similar (— 20 — 30%) decreases in collapse duration are 
seen in the 0.3 and 3 Mq initial mass cores. 

Defining the dense core star formation efficiency, fsFE, 
as the fraction of material initially in the core that ends 
up in the star at the end of collapse, this model has 
fsFE = 0.48, 0.33, and 0.31 for the 0.3, 1, and 3 Mq ini- 
tial mass cores. The higher fsFE for the 0.3 Mq initial 
mass core comes about because of the higher fraction of 
total collapse time spent in the FHSC phase (— 40%, 
compared to 5 — 10% for the 1 and 3 Mq initial mass 
cores) where there is no mass-loss. These values oi fspE 
are in general agreement with those found by Alves et 
al. (2007; 30%) by comparing the CMF of dense cores 
in the Pipe Nebula to the stellar IMF and by Enoch 
et al. (2008; > 25%) by comparing the CMF of dense 
cores in Perseus, Ophiuchus, and Serpens to the stellar 
IMF. Although the modifications described here repre- 
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sent a highly ideahzed model with representative values 
assumed for many parameters, we are encouraged by this 
agreement. 

Figure [19] compares the model 4 SEDs at various in- 
clinations for the 1 Mq initial mass core to the model 1 
SEDs presented in §4.11 There are no model 4 SEDs at 
late times since collapse ends at 165,000 yr rather than 
224,000 yr, as dicussed above. The most striking change 
compared to models 1 — 3 is the increased inclination 
dependence. Once 9c exceeds the inclination of a given 
line-of-sight the emission from the protostar and disk are 
directly observed, along with the long-wavelength emis- 
sion from the envelope. A line-of-sight 's transition from 
passing through the envelope to passing through the out- 
flow cavity is not immediate; there is a short transition 
as 9c approaches that line-of-sight 's inclination where it 
passes through both the cavity and the envelope, a result 
of the stream-line, funnel-shaped outflow cavity. An ex- 
ample of this is seen in the 58,000 and 66,000 yr panels; 
the i = 45° line-of-sight clearly shows an increase relative 
to the i = 75° line-of-sight due to geometry even though 
9c does not reach 45° until 71,000 yr. Finally, we also 
note that there is less long- wavelength emission at a given 
time compared to model 1 and that this discrepancy in- 
creases with time. Since emission at these wavelengths 
directly traces total mass, this discrepancy is due to the 
faster decrease in M^nv induced by the entrainment of 
envelope material in the outflow. 
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Fig. 19. — Same as Figure 1141 except now for model 4: SEDs 
at select ages through the collapse of the 1 Mq core. In each 
panel, the lines show the model 4 (scattering, disk, rotationally 
flattened envelope, and mass-loss and outflow cavities included) 
SEDs at various inclinations, with the line weight key given in the 
flrst panel. The dotted line shows, for comparison, the SED at that 
time for model 1. There are no model 4 SEDs at late times since 
collapse ends at 165,000 yr in the model 4 1 Mq initial mass core. 

Figures [201 [HI and [H| show the observational signa- 
tures Lboi, Tboi, and Lioi/ Lsmm plotted against the ratio 
of Mint/Mtot for the 0.3, 1, and 3 M0 initial mass cores, 
respectively. Unlike for previous models we do not com- 
bine the three masses on one plot since the increased 
inclination dependence of model 4 compared to models 
1 — 3 would create an overly complicated, difflcult-to-read 
flgure. The results discussed in relation to Figure [T9l are 
readily apparent: a given line-of-sight features low values 



of Tjjoi and Lf^oi / Lsmm until Oc approaches the inclination 
of that line-of-sight, and after a small transition region 
where both quantities increase as the line-of-sight passes 
through more of the cavity and less of the envelope, they 
increase to high values characteristic of those expected 
for viewing a protostar -|-disk directly through the outflow 
cavity. The calculated Li^d also shows an inclination de- 
pendence and is generally lower than in previous models 
due to the lower protostellar masses and lower mass ac- 
cretion rates as a consequence of including mass-loss. 




Fig. 20. — Same as Figure [TSl except now for model 4: obser- 
vational signatures Ljoj (top panels), T^oj (middle panels) and 
Lbol / Lsmm (bottom panels) versus Mint/Mtot, the ratio of in- 
ternal (protostar+disk) to total (protostar -|-disk-|-envelope) mass. 
The left panels show the model 1 results while the right panels 
show the model 4 results. The different lines show the results for 
different inclinations; the line weight key is given in the upper right 
panel. Only the 0.3 Mq model 4 core is shown to avoid creating an 
overly complicated figure. The class boundaries in T^o; are taken 
from Chen et al. (1995) while the class divisions in ijoj/L 

smm are 

taken from YE05. The discontinuities in both T^qj and itoj/Lsmm 
for model 1 are artifacts introduced by the switch from the Shu 
(1977) density profile to a power-law density profile. 

Unlike with previous models, it no longer makes 
sense to quote an inclination-weighted average value of 
Mint/Mtot when each initial mass core crosses the Class 
O/I boundary in either Thoi or Lhoi / Lsmm- Indeed, the 0.3 
M0 initial mass core crosses the Ti,oi Class /I boundary 
at values oi Mint/Mtot ranging from 0.31 to 0.91 depend- 
ing on inclination, and it crosses the Lhoi/ Lsmm Class 
O/I boundary at values of Mint/Mtot ranging from 0.33 
to never, "'^^ depending on inclination. Similar results are 
found for the 1 and 3 Mq initial mass cores. The en- 
tire concept of a connection between physical Stage and 
observational Class as measured by Tboi or Lhoi/Lsmm 
breaks down once the inclination dependence from out- 
flow cavities are taken into account, since either quantity 
can vary by an order of magnitude or more depending on 
inclination. These results are in general agreement with 
Crapsi et al. (2008), who found that Ti,oi can vary by 
factors of ~ 2 — 6 depending on the exact model pa- 
rameters. However, as their models held 9c fixed at 15° 
but adopted i — 25° as their minimum inclination, their 
maximum variation in T^oi docs not include a line-of- 
sight looking directly down the outflow cavity. With such 
small cavities, they concluded that T^oi still provided a 

^* Even when the envelope has fully dissipated, the disk keeps 
the nearly edge-on SEDs from crossing the Class O/I boundary in 
Lhol/Lsmm- The reason why this did not occur for models 2 and 
3, which also included a disk in the radiative transfer, is because of 
the lower L^oj in model 4 due to lower protostellar masses and lower 
accretion rates. Lsmm stays about the same, but L^,oi decreases, 
lowering L^^^ilLsmm in model 4 compared to models 2 and 3. 
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Fig. 21. — Same as Figure [20l except now showing the 1 Mq 
rather than the 0.3 Mq model 4 core: observational signatures Li,^i 
(top panels), T^qj (middle panels) and Lf^^i / L smm (bottom panels) 
versus Mi„t/Mtot, the ratio of internal (protostar+disk) to total 
(protostar+disk+envelope) mass. The left panels show the model 
1 results while the right panels show the model 4 results. The 
different lines show the results for different inclinations; the line 
weight key is given in the upper right panel. The class boundaries 
in Tboi are taken from Chen et al. (1995) while the class divisions 
in Lfiai/Lsmm are taken from YE05. The discontinuities in both 
Tfioi and Lt,ol/ Lsmm for model 1 are artifacts introduced by the 
switch from the Shu (1977) density profile to a power-law density 
profile. 
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Fig. 22. — Same as Figure 1211 except now showing the 3 Mq 
rather than the 1 Mq model 4 core: observational signatures L^oj 
(top panels), Tf,oi (middle panels) and Ltoi / Lsmm (bottom panels) 
versus Mi„t/Mtot, the ratio of internal (protostar+disk) to total 
(protostar+disk+envelope) mass. The left panels show the model 
1 results while the right panels show the model 4 results. The 
different lines show the results for different inclinations; the line 
weight key is given in the upper right panel. The class boundaries 
in Tjoj are taken from Chen et al. (1995) while the class divisions 
in Ljjoi/ Lsmm are taken from YE05. The discontinuities in both 
Tfjoj and Li,oi/ Lsmm for model 1 are artifacts introduced by the 
switch from the Shu (1977) density profile to a power-law density 
profile. 

good measure of physical Stage for moderate inclinations 
(25 — 70°) with lines-of-sight that do not pass through 
either the cavity or the disk. On the contrary, we find 
that neither T^oi nor L^oi / Lsmm provides a good measure 
of physical Stage regardless of inclination. 

Figure [221 shows a BLT diagram for model 4. Unlike 
with previous models, the full extent of the embedded 
sources in L^oi — T^oi space is reproduced by model 4. 
However, Figure [Ml which plots model 4 Ltoi and Tt,oi 
histograms, shows that while the distribution of time 
spent at various luminosities is wider in model 4 than 
in models 1 — 3 and clearly gives a better fit to the ob- 
served distribution (a K-S test gives a 22% probability 
that the observed and model Lboi histograms represent 
the same underlying distribution, compared to < 0.1% 
for models 1 — 3), the model still overpredicts the time 
spent at ~ 2 — 20 Lq and underpredicts the time spent at 
'--^ 0.1 — 2 Lq. Figure [23l also shows that the model spends 



a relatively large fraction of time at high Tboi 1000 K) 
compared to the fraction of embedded sources at such 
values, a consequence of viewing direct protostar+disk 
emission though outflow cavities for many model lines- 
of-sight. This is also evident in both Figure [24l (a K-S 
test gives a 16% probability that the observed and model 
Tboi histograms represent the same underlying distribu- 
tion) and column 7 of Table [H which lists the fraction 
of total time model 4 spends in various L{,oi — T^oi bins. 
The model spends 40.1% of the time at Tboi > 1000 K 
whereas only 4.5% of the embedded sources are found at 
such high Tboi- The model spends most of the rest of the 
time at low Tboi 100 K) whereas the embedded obser- 
vations are relatively evenly distributed (in a logarithmic 
binning) between '--^ 20 — 500 K. 
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Fig. 23. — Same as Figure 1161 except now for model 4. The 
grayscale pixels indicate the fraction of total tim e the model spends 
in each Ljoj — T;,;,; bin, calculated from Equation llOl The grayscale 
is displayed in a logarithmic stretch with the scaling chosen to 
emphasize the full extent of the models in Lf^^i — Ti,^i space. The 
mapping between grayscale and fraction of total time is indicated 
in the legend. The class boundaries in Tijoj are taken from Chen 
et al. (1995). The thick dashed line shows the ZAMS (D'Antona 
& MazziteUi 1994) from 0.1 to 2.0 Mq. The colored symbols show 
the Young Stellar Objects from Evans et al. (2009) in this diagram; 
the colors and symbols hold the same meaning as in Figure [7] 
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Fig. 24. — Same as Figure lTTl except now for model 4: histograms 
showing the fraction of total sources (observations; solid lines) and 
fraction o f to tal time spent (models; dashed lines; calculated from 
Equation IIOI I at various L^oj (left) and Tjoj (right). The binsize 
is 1/3 dex in both quantities. For the observations, only the 112 
embedded sources (plotted as filled circles on the BLT diagrams) 
are included. 

In summary, the main conclusions of models 1 — 3 that 
the model overpredicts both the time spent at high Lboi 
(> 1 - 2 Lq) and the time spent at low Tboi (< 100-200 
K) remain unchanged. However, including the effects of 
mass-loss and outflow cavities has reduced the severity of 
the luminosity problem and also introduced a significant 
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population of model 4 SEDs with higher Ttoi than found 
for embedded sources. 

Before moving on, we briefly return to the assumed 
value of /, the efficiency with which the jet /wind trans- 
fers its momentum to the ambient medium. We assumed 
representative values based on either observed or theo- 
retical ranges for all other relevant parameters, but we 
maximized / by setting it equal to 1. We made this 
choice to maximize the effects of mass-loss and outflow 
cavities, since / = 1 maximizes the amount of mass en- 
trained in the outflow and thus the increase in 0c with 
time. Even with their effects maximized, mass-loss and 
outflow-cavities still feature the same basic luminosity 
problem as the other models, albeit to a lesser degree. 
We consider this to be a strong test of the necessity of 
invoking episodic accretion to explain the observed distri- 
bution of embedded sources in Li,oi — T^oi space. How- 
ever, with f = 1, 6c reaches 90° before collapse ends 
(defined as tcoUapse = M*™**'^7M*™**"') and thus termi- 
nates the embedded phase earlier than when no outflow 
cavities are present, and the large cavities prevent both 
Tfioi and Lhoi/ Lgmm from being useful indicators of physi- 
cal Stage. What if we had assumed a more common value 
in the range of / = 0.1 - 0.25 (e.g., Andre et al. 1999)? 

With / = 0.1, collapse ends before 9c reaches 90° in all 
three initial mass cores. However, these angles still reach 
60, 65, and 75° at the end of the collapse of the 0.3, 1, 
and 3 Mq initial mass cores, respectively. Thus, we still 
disagree with the conclusion of Crapsi et al. (2008) that 
Tboi provides a good measure of evolutionary status for 
moderate inclinations (25 — 70°) with lines-of-sight that 
do not pass through either the cavity or the disk, which 
they reached by assuming a constant dc — 15°. Fur- 
thermore, the measured values of fsFE with / = 0.1 are 
0.72, 0.64, and 0.64, respectively, much higher than those 
measured with / = 1 and found above to be in general 
agreement with observational estimates of this quantity. 
Finally, / = 0.1 means that the protostellar mass will 
grow more quickly and the envelope accretion rate will 
decrease more slowly, both of which will increase model 
luminosities that are either dominated by or include sig- 
nificant components from accretion luminosity (see ji2.2[) . 
We thus conclude that choosing more common values of 
0.1 — 0.25 for the entrainment efficiency will not change 
any of our basic results, and will lessen the degree to 
which the luminosity problem is improved by model 4. 

4.5. Model 5: Episodic Accretion 

In models 1 — 4 we have modified the YE05 model 
to include isotropic scattering off dust grains, a two- 
dimensional disk in the radiative transfer, rotationally 
flattened envelope density profiles following the TSC84 
solution for the collapse of slowly rotating cores, and 
mass-loss and outflow cavities. All of these have had 
impacts on the model predictions, especially including 
the opacity from scattering and mass-loss and outflow 
cavities. However, even with all the above effects consid- 
ered, the models still exhibit the fundamental luminos- 
ity problem described in SJTJ the models overpredict the 
fraction of total time spent at high luminosities (> 1 — 2 
Lq) compared to observations of embedded protostars 
(although the effects of mass-loss and outflow cavities re- 
duce the severity of the problem). Given this, along with 



the discussion in fJT] regarding observational evidence for 
non-steady mass accretion and theoretical predictions of 
episodic mass accretion, we incorporate episodic mass 
accretion into our evolutionary model. As with ^ 14.41 we 
incorporate a simple, idealized process of episodic accre- 
tion that, while physically motivated, is designed to test 
its general effects on the evolutionary signatures of em- 
bedded protostars rather than fully capture a complete, 
physically self-consistent model. 

The physical basis for our treatment of episodic ac- 
cretion are the models published by Vorobyov & Basu 
(2005b, 2006), who use MHD simulations to follow the 
collapse of rotating cores. In their simulations, material 
piles up in a circumstellar disk until the disk becomes 
gravitationally unstable and develops spiral structure 
and dense clumps, which are then driven onto the pro- 
tostar in short-lived accretion bursts generated through 
the gravitational torques associated with the spiral arms. 
They found this burst phenomenon to be a robust result 
under a variety of initial conditions (Vorobyov & Basu 
2006). Other authors have also shown that gravitational 
instabilities in the disk can produce episodic mass accre- 
tion onto the protostar. For example, Boss (2002) noted 
that their models of gravitationally unstable disks (used 
primarily to investigate giant planet formation) feature 
protostellar mass accretion rates that vary with time be- 
tween about 10"'^ to 10"^ Mq yr'^ 

To incorporate episodic accretion into our model, we 
allow the envelope to evolve as before: material accretes 
from the envelope directly onto the protostar and disk 
(but mostly the disk as it grows in size), with the enve- 
lope density profile evolving following the TSC84 solu- 
tion for the collapse of slowly rotating cores. This leaves 
the 1st, 2nd, and 4th luminosity sources as described 
in unchanged. However, instead of material accret- 
ing through the disk and onto the protostar at the rate 
MotoP — VoMcnv, all material accreted onto the disk is 
stored in the disk {MotoP is set to zero). This causes the 
3rd luminosity source, that arising from accretion from 
the disk onto the protostar, to vanish, and the 5th lu- 
minosity source, which depends on the total mass accre- 
tion rate onto the protostar, to be significantly reduced. 
Since the luminosity arising from accretion from the disk 
onto the protostar is the dominant source of luminosity 
throughout most of the model evolution (due to the deep 
potential well of the protostar), this will significantly re- 
duce the total model luminosity. 

The disk mass can't continue to grow indefinitely; 
eventually gravitational instabilities will set in once its 
mass grows to a significant fraction of the protostellar 
mass. In a study of disks around Class II sources in Ophi- 
uchus and Taurus, Andrews & Williams (2007) found 
a Md/M^ distribution ranging from ^ 10^"^ — 0.2. If 
we assume the upper end of this distribution represents 
the youngest disks that have not yet begun to disperse, 
this should be a good estimate of the maximum possible 
Md/M^ as any disk higher than this would have become 
gravitationally unstable. This is in good agreement with 
theoretical predictions of when gravitational instabilities 
develop (Md/M, 0.2; e.g., Shu et al. 1990), and also 
with Vorobyov & Basu (2006), who showed that the disk 
mass always remains significantly less than the protostel- 
lar mass in their simulations. More recently, Vorobyov 
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(2009a) argued that the maximum Md/M^, is closer to 
0.25 — 0.4; such a higher value would increase the time 
between bursts but also the duration of each burst, and 
would not significantly change our results. 

Thus, once M^/M^, reaches 0.2, the model enters a 
burst mode and MotoP increases to 1 x 10"'* Mq yr~^. 
In their simulations, Vorobyov & Basu (2005b, 2006) 
found that MotoP increases to about (1 — 5) xl0~^ Mq 
yr""'^, even reaching ~ 10~^ yr~* in the most ex- 
treme cases; here we assume a constant value of 1 x lO"'' 
Mq yr"""^ for every burst for simplicity. In the bursts we 
increase the time resolution from A< = 1000, 2000, or 
6000 yr for the 0.3, 1, and 3 Mq initial mass cores, re- 
spectively, to At = 100 yr. The high MotoP cause s bo th 
the 3rd and 5th luminosity sources as described in §2.2l to 
increase to very high values, dominating the total lumi- 
nosity during bursts (see below). The protostellar mass 
grows rapidly during a burst, and the disk mass decreases 
accordingly (during a burst the material is accreting out 
of the disk about two orders of magnitude faster than it 
is accreting onto the disk from the envelope). For sim- 
plicity the burst is assumed to continue until all of the 
disk mass accretes onto the protostar, at which point 
MotoP drops back to zero, the time resolution decreases 
back to its original value, and the cycle begins anew. We 
assume the dust temperature has reached equilibrium by 
the timesteps immediately following the onset and ter- 
mination of accretion bursts*^. 

Figure [25] shows the effects of including episodic accre- 
tion as described above for the 1 Mq initial mass core. 
Similar to Figure fTSl for model 4, plotted are the masses 
of the various model components (protostar, disk, enve- 
lope, and outflow), the cavity semi-opening angle, and 

Menv versus time. Except for the first 0.02 Myr when 
the disk has not yet formed, shows a "staircase" func- 
tion, essentially increasing only during the bursts (except 
for the very small amount of accretion directly from the 
envelope onto the star). Each increase in is accompa- 
nied by a corresponding decrease in M^. Each burst also 
features an increase in the outflow mass and decrease in 
Menv since mass is ejected (and thus envelope mass is en- 
trained) when material accretes from the protostar to the 
disk. These increases in the outflow mass are accounted 
for by increases in the outflow cavity semi-opening angle 
and thus cause decreases in Menv following Equation [TBI 
The collapse ends at 180,000 yr; a burst begins at this 
time and by the next timestep, 100 yr later, the outflow 
cavity semi-opening angle has reached 90°, removing the 
remaining envelope material and ending collapse. 

The 0.3, 1, and 3 Mq initial mass cores feature 3, 8, and 
11 bursts, respectively. Vorobyov & Basu (2005b) found 
15 — 30 bursts in their initial simulation, but showed 
in Vorobyov & Basu (2006) that the exact number de- 
pends on assumed values of both and the magnetic 

Once the change in luminosity reaches the dust, it will respond 
very quickly to the change, reaching its new equilibrium within 
a few seconds (e.g., Draine &: Anderson 1985; Lee et al. 2007). 
The random-walk time through the envelope is the hmiting factor, 
and a conservative upper hmit calculated assuming photons remain 
at their initial wavelength instead of being reprocessed to longer 
wavelengths (where the optical depth is lower) puts this timescale 
at ~ 100 yr. In reality it is much less once reprocessing to longer 
wavelengths is considered. 




Time (Myrl 

Fig. 25. — Same as Figure [TSl except now for model 5: evolution 
of the masses of various model components (top panel), outflow 
cavity semi-opening angle (middle panel), and envelope mass ac- 
cretion rate (bottom panel) versus time for the 1 NLq initial mass 
model 5 core, which includes episodic accretion as described in i|4.5l 

field. It will also depend on the exact criterion assumed 
for gravitational instability, which we have held fixed at 
Mrf/M* = 0.2 for simplicity. All three initial mass cores 
spend 1.5 — 2% of their total collapse times in burst 
phases and have fsFE of 0.53, 0.35, and 0.33, respec- 
tively. As expected, these are similar to model 4 since 
including episodic accretion changes when material ac- 
cretes onto the protostar but does not generally affect 
the total amount of material accreted. 

Showing SEDs at select ages through the collapse of 
the 1 Mq initial mass core is less meaningful here than 
for models 1 — 4 since the evolution from starless core to 
revealed protostar-|-disk is no longer smooth but changes 
abruptly in the bursts. We instead show, in Figure (261 
the SEDs just before and during accretion bursts that 
begin at t = 24400 yr and t = 100800 yr. While the core 
itself is at a much different evolutionary state at the two 
times, the effects of an accretion burst are similar: the 
flux increases due to the increase in luminosity, and the 
overall shape of the SED shifts to shorter wavelengths 
due to the increase in emission from both the protostar 
itself and the warm dust in the envelope heated by the 
protostar. 

Figures [27l [28l and [29] show the observational signa- 
tures Lboh Tboh and Lboi/ Lsmm plotted against the ratio 
of Mint/Mtot for the 0.3, 1, and 3 Mq initial mass cores, 
respectively. As with model 4, we show each initial mass 
core in a separate figure to avoid creating a confusing, 
difficult-to-read figure. Since the envelope mass quickly 
decreases in each burst due to the entrainment of ma- 
terial in the outflow (see above and Figure [25]), Mtot 
also quickly decreases in each burst and thus the ratio of 
Mint/Mtot quickly increases. As a result, the values on 
the X-axis in Figures [27] — [22] no longer increase linearly 
with time; the bursts occupy a much lower fraction of 
total time than they do total Mint/Mtot- 

Without accretion from the disk onto the protostar 
providing the dominant source of luminosity, most of the 
time is spent at low luminosity (often > 1 — 2 orders of 
magnitude lower compared to model 1) except during the 
bursts, when the luminosity increases to very high values 
100 — 1000 Lq). The general conclusion from model 
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Fig. 26. — Select SEDs in the collapse of the model 5 IMq core. 
The top panels show the SEDs just before (t = 24000 yr; right) 
and during {t = 24400 yr; left) an accretion burst that begins 
at t = 24400 yr. The bottom panels show the SEDs just before 
(t = 100000 yr; right) and during {t = 100800 yr; left) an accretion 
burst that begins at i = 100800 yr. Three different inclinations are 
shown: nearly pole-on (i = 15°), moderate (i = 45°), and nearly 
edge-on (i = 75°), with the line weight key given in the first paneh 
In the top two panels all three inclinations feature nearly equiv- 
alent SEDs since the departure from spherical symmetry is very 
small at early times. In the bottom two panels the nearly pole-on 
and moderate inclinations feature nearly equivalent SEDs because 
they are both viewing direct protostar-pdisk emission through the 
outflow cavity. 




Fig. 27. — Same as Figure 1201 except now for model 5: obser- 
vational signatures L^qj (top panels), T(,oj (middle panels) and 
Lbol / Lsmm (bottom panels) versus Mint/Mtot, the ratio of in- 
ternal (protostar+disk) to total (protostar-|-disk-|-envelope) mass. 
The left panels show the model 1 results while the right panels 
show the model 5 results. The different lines show the results for 
different inclinations; the line weight key is given in the upper right 
panel. Only the 0.3 Mq model 5 core is shown to avoid creating an 
overly complicated figure. The class boundaries in Tj,oj are taken 
from Chen et al. (1995) while the class divisions in Li,ol/L smm are 
taken from YE05. The discontinuities in both Tjoj and Li,ai/Lsmm 
for model 1 are artifacts introduced by the switch from the Shu 
(1977) density profile to a power-law density profile. 



4 that the large inchnation dependence introduced by 
the outflow cavities prevents both T^oi and L{,oi/ Lgmm 
from providing good measures of evolutionary status re- 
mains unchanged. For example, excluding time spent 



in bursts, the 1 Mq initial mass core crosses the T; 
Class /I boundary at values of Mint /Mtot ranging from 
0.11 — 0.90 depending on inclination. Furthermore, again 
excluding time spent in bursts, this core spends essen- 
tially all of its time as Class by L^oil Lsmm, crossing 
the Class O/I boundary at values of Mint /Mtot ranging 
from 0.84— never depending on inclination. Both quan- 
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Fig. 28. — Same as Figure 1281 except now showing the 1 Mq 
rather than the 0.3 Mq model 5 core: observational signatures Lf,„i 
(top panels), T^oj (middle panels) and Lf^^i / L smm (bottom panels) 
versus Mi„t/Mtot, the ratio of internal (protostar-|-disk) to total 
(protostar+disk-|-envelope) mass. The left panels show the model 
1 results while the right panels show the model 5 results. The 
different lines show the results for different inclinations; the line 
weight key is given in the upper right panel. The class boundaries 
in Tfjoj are taken from Chen et al. (1995) while the class divisions 
in Li,ai/Lsmm are taken from YE05. The discontinuities in both 
Tijoj and Li,ol / Lsmm for model 1 are artifacts introduced by the 
switch from the Shu (1977) density profile to a power-law density 
profile. 




Fig. 29. — Same as Figure 1281 except now showing the 3 Mq 
rather than the 1 Mq model 5 core: observational signatures L^oj 
(top panels), Ttoi (middle panels) and Ltoi / Lsmm (bottom panels) 
versus Mi„t/Mtot, the ratio of internal (protostar-|-disk) to total 
(protostar+disk-|-envelope) mass. The left panels show the model 
1 results while the right panels show the model 5 results. The 
different lines show the results for different inclinations; the line 
weight key is given in the upper right panel. The class boundaries 
in Tfjoj are taken from Chen et al. (1995) while the class divisions 
in Liiol/ Lsmm are taken from YE05. The discontinuities in both 
Ttot and Li,ol/Lsm?n for model 1 are artifacts introduced by the 
switch from the Shu (1977) density profile to a power-law density 
profile. 

titles are similarly unreliable at measuring the evolution- 
ary status of the 0.3 and 3 Mq initial mass core. 

However, episodic accretion adds another layer of unre- 
liability to these evolutionary indicators. Both quantities 
clearly show large increases during bursts and subsequent 
decreases after each burst ends. Depending on the com- 
bination of inclination and age, Thoi and Li,oi/Lsmm cal- 
culated for a given line-of-sight can cross back and forth 
across the Class O/I boundary in both quantities several 
times throughout the collapse of the core. If this model 
represents physical reality, which we evaluate below, nei- 
ther of the commonly used evolutionary indicators actu- 
ally reliably trace physical Stage. 

Figure [501 shows a BLT diagram for model 5. The full 
extent of the embedded sources in Lboi — Tboi space is 
again reproduced, as for model 4. The model reaches 
higher luminosities than models 1 — 4 as a result of the 
accretion bursts. At first glance these high luminosi- 
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ties appear inconsistent with observations, although it is 
difficult to evaluate from this figure alone given the log- 
arithmic grayscale. The model only spends ~1% of its 
total time at Lboi > 100 L©;^^ for comparison, none of 
the observed sources have L^oi > 100 L©. Given that the 
observed sample only consists of 112 sources, ^--^1% only 
corresponds to about one source, and thus it is not in- 
consistent, given the small-number statistics, to find no 
sources at such high luminosities. 
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Fig. 30. — Same as Figure 1231 except now for model 5. The 
grayscale pixels indicate the fraction of total tim e the model spends 
in each Ljoj — Tjoj bin, calculated from Eguation llOl The grayscale 
is displayed in a logarithmic stretch with the scaling chosen to 
emphasize the full extent of the models in Lt,ol — Ttot space. The 
mapping between grayscale and fraction of total time is indicated 
in the legend. The class boundaries in Tjoj are taken from Chen 
et al. (1995). The thick dashed Une shows the ZAMS (D'Antona 
& MazziteUi 1994) from 0.1 to 2.0 Mq. The colored symbols show 
the Young Stellar Objects from Evans et al. (2009) in this diagram; 
the colors and symbols hold the same meaning as in Figure [3 

Figure [nH plots model 5 L^oi and Tboi histograms, and 
column 8 of Tablc[l]lists the fraction of total time model 5 
spends in various Ltoi — Tboi bins. A K-S test shows that 
there is a 61% probability that the observed and model 
Lboi histograms represent the same underlying distribu- 
tion, by far the highest of models 1 — 5. Inspection of 
Figure [3T] shows that most of the remaining discrepancy 
between the models and observations is in the form of a 
sort of "reverse luminosity problem" in that the models 
spend an excess of time at Lboi ~ 0.1 Lq compared to 
the observations. In this model we assumed a mass ac- 
cretion rate from the disk to the protostar of MntoP = 
Mq yr~^ in between accretion bursts. In reality, how- 
ever, accretion from the disk onto the protostar likely 
continues at a low rate during the quiescent phases; in 
their simulations Vorobyov & Basu (2005b, 2006) found 
that MotoP is typically 10"^ - 10"^ M© yr^^ dur- 
ing these phases. Such non-zero values of MutoP, even 
though they are low, could easily add a few xO.l L© to 
the total luminosity since this term dominates the to- 
tal luminosity when present (see ji5.4[) . Thus, we do not 
consider the excess of time spent at Lboi ~ 0.1 L© in the 
model to be significant, and we conclude that model 5 
best reproduces the observed luminosity distribution of 
embedded protostars and is the only model that essen- 
tially resolves the luminosity problem. 

^® Even though each core spends about 2% of the total collapse 
time in burst, and the total model luminosity in every burst is > 
100 L0, the higher inclination lines-of-sight that are weighted more 
heav ily due to the increased solid angle at higher inclinations [see 
i|4.1| often have Lh^i < 100 Lq even in accretion bursts 
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Fig. 31. — Same as Figure [24l except now for model 5: histograms 
showing the fraction of total sources (observations; solid lines) and 
fraction o f to tal time spent (models; dashed lines; calculated from 
Equation I10| l at various L),oi (left) and T(,oi (right). The binsize 
is 1/3 dex in both quantities. For the observations, only the 112 
embedded sources (plotted as filled circles on the BLT diagrams) 
are included. 

Finally, we note that a K-S test gives a 29% probability 
that the observed and model Thoi histograms represent 
the same underlying distribution, comparable to models 
1 — 4. As with model 4, model 5 includes a significant 
population of SEDs with higher Thoi than found for em- 
bedded sources (35.6% of the time is spent at Tbd > 1000 
K whereas only 4.5% of the embedded sources are found 
at such high Tjyoi)- Again, this is a consequence of the 
outflow cavities allowing many lines-of-sight to view di- 
rect protostar-|-disk emission, increasing the 1 — 100 /im 
emission and thus the calculated Tjjoi- We will return to 
this point in Section [5.21 

5. DISCUSSION 
5.1. Scattering and 2-D Geometry 

Including the opacity from scattering has an important 
effect on the results. Since the opacity from scattering 
dominates over that from absorption at approximately 
the same wavelengths over which the protostar-|-disk in- 
put SED peaks (~ 0.1 — 10 ^m), including this opacity 
significantly increases the total optical depth through the 
model, causing a corresponding decrease in the amount of 
emission at these wavelengths escaping from the model. 
This reduction in the short-wavelength emission causes 
a decrease in measured values of Tboj (by factors of 
~ 1.5 — 6) and slows down the increase in T^oi as the 
model evolves. As a consequence, the model crosses the 
Class O/I boundary in Tjjoi at later times than found by 
YE05. We disagree with their conclusion that T^oi is not 
a good evolutionary indicator; we find that it is not quite 
as good as L^oi / Lsmm but remains a satisfactory indica- 
tor of evolutionary status (although this will change in 
the later models; see below). Despite the important ef- 
fects of including the opacity from scattering, the general 
luminosity problem described in fJT] remains. 

Switching to a 2-D model setup and including a cir- 
cumstellar disk and rotationally flattened envelope den- 
sity profile in the model both introduce an inclination de- 
pendence, as expected. However, the dependence is only 
significant at late times in the collapse of the cores (both 
because the disk only becomes relatively large and mas- 
sive at late times, and because the degree of flattening 
of the envelope starts out very small and increases with 
time). Even at such late times the inclination depen- 
dence, while introducing spread in the evolution of Lboi, 
Tboi, and Lboi/ Lsmm with time, is not large enough to 
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significantly change any of the conclusions (see ^4.21 and 
14.31 for a quantitative description of the changes). Thus, 
while important and more physically realistic, these ge- 
ometry effects do not resolve the discrepancy between 
the model and the observations. 

5.2. Mass-Loss, Outflow Cavities, and Episodic 
Accretion: Connection Between Observational 
Class and Physical Stage 

Including the effects of mass-loss, outflow cavities, and 
episodic accretion has significant effects on the observa- 
tional signatures of the model. As shown in both ^4:A\ 
and !j4.51 the variation in T^oi and L^oi / Lsmm with incli- 
nation is dramatically increased by the outflow cavities 
as lines-of-sight with inclinations less than 6c (the cavity 
semi-opening angle) view direct protostar-|-disk emission 
through the cavity in addition to the long-wavelength 
emission from the envelope. T^oi and L{,oi/ Lsmm calcu- 
lated from SEDs at the same time in the collapse of the 
same initial mass core can vary by an order of magni- 
tude or more depending on inclination, and both quan- 
tities cross their Class O/I boundaries at times ranging 
from very early in the collapse for nearly pole-on lines 
of sight to very late in the collapse (or even sometimes 
never for L^oi/ Lgmm) for nearly edge-on lines of sight. 
This led us to conclude in both ^14. 41 and that obser- 
vational Class determined by Thoi or Lj^oi/ Lsmm is not a 
good indicator of physical Stage. To quantify this. Table 
[3] lists the fraction of total time that model 5 spends in 
the various Class/Stage combinations according to both 
Tboi and Lboi/ Lgmm- From this table it is apparent that 
Class determined by Ttoi only agrees with the underlying 
physical Stage 40% of the time; the remaining 60% is oc- 
cupied by either Ti,oi measuring Class II while still in the 
embedded phase, or Tboi measuring Class or I while the 
true physical Stage is the opposite. L^oi/ Lgmm, on the 
other hand, gives a Class that agrees with the underly- 
ing physical Stage 79% of the time. Class determined by 
Lboi I Lsmm is thus more likely to trace the true physical 
stage of an embedded object than if it is determined by 
Tboi, although we caution that the agreement would be 
worse if Lboi/Lsmm included a defined Class I/II bound- 
ary. We also caution that Lboi/Lsmm is more difficult to 
determine accurately since the calculated value is signifi- 
cantly more sensitive than T^oi to the exact wavelengths 
at which one samples the SED with observations (Dun- 
ham et al. 2008). Thus, in general, neither evolutionary 
indicator provides a reliable measure of physical Stage, 
although Lboi/ Lsmm is statistically more likely to pro- 
vide an accurate measure of physical Stage than Ttoi- 
Furthermore, the entire concept ofTboi or L^oilLsmm in- 
creasing monotonically and progressing through the var- 
ious Classes as the object evolves from Stage through 
Stage I to Stage II at the end of the embedded phase 
breaks down once episodic accretion is included. 

This point was recognized by Myers et al. (1998), 
who acknowledged that the effects of geometry can limit 
the utility of T^oi- However, the variation in T^oi and 
Lboi/Lsmm with inclination is generally larger in this 
model than found by previous authors, mainly because 
we allow 9c to increase to very large values (eventually 
reaching 90° and ending collapse). For example, Crapsi 
et al. (2008) held 9c fixed at 15° but only presented 
results for i > 25°, effectively eliminating any lines-of- 



sight looking directly down the outflow cavity. Myers 
et al. (1998) concluded that outflow cavities will only 
change Tboi by about a factor of 2 for 0c = — 25°, 
which they assumed to be representative of typical cavi- 
ties, although they did acknowledge the variation will be 
larger for larger cavities. Whitney et al. (2003a, 2003b) 
assumed small outflow cavities {9c ~ 20 — 25°) while 
Robitaille et al. (2006) allowed cavities as large as 60°; 
both showed that classification according to infrared col- 
ors can vary with inclination but neither discussed the 
effects on Tboi or Lboi/Lgmm- 

Are the outflow cavities featured here, with 9c increas- 
ing to such large values, consistent with observations? 
Based on radiative transfer models, Tobin et al. (2007) 
found 9c = 15 — 20° for 3 Class sources, and Furlan 
et al. (2008) found 9c = 0.1 - 27° for 22 Class I sources. 
However, all 25 sources in the combined sample have 
9c < i, and as Furlan et al. point out, the lack of any 
sources with large 9c could be a selection bias introduced 
by the fact that such sources may be classified as Class 
II and thus missing from the sample. Indeed, some ob- 
servational evidence for larger outflow cavities is found 
by Huard et al. (2006), who measured 9c ~ 50° for the 
Class very low luminosity object (VeLLO; Di Francesco 
et al. 2007) L1014 from deep near-infrared imaging, and 
by Arce & Sargent (2006), who measured 6'c ~ 10 — 60° 
for a sample of 6 Class O/I embedded sources from mil- 
limeter interferometer spectroscopy. 

Thus, there is some evidence that outflow cavities can 
be larger than assumed by Crapsi et al. (2008) and Myers 
et al. (1998). By choosing an entrainment efficiency be- 
tween the ejected jet/wind and ambient medium of 100%, 
we maximized the speed with which 9c increases in order 
to maximize the effects of mass-loss and outflow cavities 
and test whether or not episodic accretion is truly needed 
to resolve the luminosity problem (see below). However, 
as discussed in tj4.41 even a much smaller, more typi- 
cally assumed value of 10% for the entrainment efficiency 
opens cavities that reach 9c = 60, 65, and 75° at the end 
of the collapse of the 0.3, 1, and 3 M© initial mass cores, 
respectively. In addition to the fact that such a model 
would lessen the degree to which the luminosity problem 
is improved (discussed in i i4.4l and more detail below) , 
it also gives star formation efficiences higher than those 
found when the assumed entrainment efficiency is 100%, 
which we showed in ^AA\ are generally consistent with es- 
timated values from observations. More fundamentally, 
there is a general inconsistency between estimated val- 
ues of the star formation efficiency (~ 30%; e.g., Alves 
et al. 2007; Enoch et al. 2008) and small outflow cavi- 
ties. Assuming the cavities reached their maximum sizes 
immediately after collapse begins, only 13 — 50% of the 
mass would be removed for 9c = 30 — 60°. Simple conser- 
vation of momentum and the velocity difference between 
the jet/wind and outflow argues that the bulk of the mass 
removed from the system is in the outflow rather than 
the jet/wind, thus such cavities don't remove enough ma- 
terial to match current estimates of the star formation 
efficiency. The assumption of the formation of more than 
one star per core could help to alleviate this discrepancy, 
but, on the other hand, outflow cavities don't reach their 
maximum sizes immediately, making the 13 — 50% an up- 
per limit only. 

The outflow cavities create a significant population of 
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embedded objects at high Tboi (> 1000 K). Is it possible 
that a population of embedded objects with such high 
Thoi exists? No such population is included in the Evans 
et al. (2009) sample of embedded objects. This sample 
is based primarily on the work of Enoch et al. (2009) 
and Dunham et al. (2008); the former compiled a com- 
plete sample of embedded objects in Perseus, Serpens, 
and Ophiuchus while the latter presented the results of 
a search for all low- luminosity (< 1 Lq) embedded ob- 
jects in the molecular clouds and isolated cores observed 
by the Spitzer Space Telescope Legacy Project "From 
Molecular Cores to Planet Forming Disks" (Evans et 
al. 2003). While both studies devoted considerable at- 
tention to completeness, many of the criteria for identi- 
fying candidate embedded objects were based on detec- 
tions at 24 — 70 fim and source SEDs exhibiting rising 
fluxes from shorter to longer wavelengths. As all panels 
at t > 34000 yr in Figure [19] clearly show, lines-of-sight 
passing through outflow cavities and seeing direct pro- 
tostar-|-disk emission often do not exhibit rising fluxes 
at 2 — 70 /j,m and do not generally show SEDs typically 
associated with embedded objects. It is possible such a 
population of embedded objects does exist but is not in- 
cluded in the Evans et al. (2009) Class I sample because 
of selection biases. It is difficult to evaluate the exact 
criteria assumed by Enoch et al. and Dunham et al. to 
determine if the fraction of time the model spends at 
high Tboi would be recovered in their samples since both 
studies make use of an automatic source classification 
scheme only available for sources observed by the Spitzer 
Space Telescope (see Evans et al. 2007 for details) and 
both included significant human judgment to determine 
what was and was not an embedded source. 

Uncertain extinction corrections may also play a role 
in this discrepancy between observations and models. 
When comparing observations of embedded sources to 
models, corrections must be applied to remove fore- 
ground extinction arising from the molecular cloud and 
ISM (separate from local extinction by the envelope, 
which will be reradiated in the far-infrared). Evans et 
al. (2009) did correct the 112 embedded sources that are 
plotted on the BLT diagrams and used to make the his- 
tograms in this paper for foreground extinction. How- 
ever, in practice it is difficult to determine the value of 
the foreground extinction to an embedded protostar, so 
Evans et al. simply applied the mean extinction to all 
the Class II objects in the same cloud. If protostars 
form in denser parts of clouds, it is possible most embed- 
ded sources are undercorrected for foreground extinction. 
To quantify this, we took the 1 Mq initial mass core at 
150,000, extincted the model SEDs by Ay = 10 and 20, 
and then un-extincted them (corrected for extinction) as- 
suming Ay = 5.9 (the mean value for Perseus, which con- 
tributes approximately half of the total sample; see Evans 
et al. [2009]). This decreases Tboi for hnes-of-sight look- 
ing through the outflow cavity from ^3000 K to ~1600 
K for Ay = 10 and ~450 K for Ay = 20. Future work 
must revisit the observational samples and carefully eval- 
uate whether or not a population of embedded sources 
viewed through outflow cavities and thus exhibiting high 
Tboi and Lboi/Lsmm exists, and also whether or not calcu- 
lated Tboi and Lboi / Lgmm suffer from an undercorrection 
for foreground extinction. 



5.3. Mass-Loss, Outflow Cavities, and Episodic 
Accretion: Towards Resolving the Luminosity 
Problem 

The primary motivation for including the modifica- 
tions to the YE05 model in the step-by-step fashion de- 
scribed in was to test the hypothesis that episodic 
accretion is necessary to resolve the luminosity problem 
and explain the distribution of sources in Lboi — Tboi 
space by eliminating other possibilities. While each had 
important effects and improved the physical realism of 
the model, including the scattering from opacity and 2- 
D effects of a circumstellar disk and rotationally flattened 
envelope left the conclusions essentially unchanged: the 
model spent too much time at high Lboi {■^ ^ ~ '2 Lq) 
and low Tboi 100 — 200 K) compared to that expected 
from the distribution of embedded sources in Lboi — Tboi 
space. 

Including the effects of outflow cavities and mass-loss 
lessened the severity of the luminosity problem but did 
not eliminate it, even when the effects were maximized 
by assuming a 100% entrainment efficiency between the 
jet/wind ejected by the protostellar system and the sur- 
rounding envelope. A smaller, more typically assumed 
value of 10% lessened the degree to which the luminosity 
problem was resolved (and also increased the star for- 
mation efficiencies to values higher than expected from 
observationally determined estimates). Thus, even with 
mass-loss maximized, which minimizes both the proto- 
stellar mass and the mass accretion rate and thus mini- 
mizes the model luminosity, the model still overpredicts 
the amount of time spent at Lbd > 1 — 2 Lq . We consider 
this to be a strong indication of the necessity of invok- 
ing episodic accretion to bring models in agreement with 
observations, on top of the observational and theoretical 
evidence for such a process described in ^ 

Indeed, model 5, which includes a simple treatment 
of episodic mass accretion based loosely on the simula- 
tions by Vorobyov & Basu (2005b, 2006) on top of the 
other modifications, is the only model that essentially 
resolves the luminosity problem. If episodic accretion 
does in fact occur, as supported by the models presented 
in this paper, there may be important consequences for 
planet formation since the properties of the circumstellar 
disk at the end of the embedded stage, in particular the 
disk mass, will depend on where in the cycle of episodic 
accretion the system is when the envelope fully dissi- 
pates. Another consequence is that the accretion bursts 
account for a large fraction of the total accretion onto 
the protostar: 50%, 83%, and 91% of the final stellar 
mass accretes during these bursts for the model 5 0.3, 1, 
and 3 Mq initial mass cores. The range is due almost 
entirely to different fractions of total time occupied by 
the FHSC phase (the first 20,000 yr) where mass accretes 
directly from the envelope onto the protostar. Thus, if 
this simple model refiects reality as comparison to obser- 
vations suggests, between 50-90% of the final protostellar 
mass accretes in < 2% of the total duration of the em- 
bedded phase. This is in general agreement with Evans 
et al. (2009), who used their luminosity distribution of 
embedded sources and a simple toy model to conclude 
that 50% of the final protostellar mass accretes in 7% of 
the lifetime. As they noted, 7% is an upper limit only 
since their sample may lack the rarest, most luminous 
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accretion events, a suggestion reinforced by our results. 

Finally, we caution that these results do not prove 
episodic accretion occurs, either as described by 
Vorobyov & Basu (2005b, 2006) or in some other fash- 
ion. While we consider the results of this paper to be 
strong evidence in favor of a process of episodic accretion 
existing in the formation of low-mass protostars, future 
work must continue to search for definitive observational 
evidence that such a process occurs (see discussion in 
Section 6.4 of Dunham et al. [2008]). 

5.4. Model Assumptions 

The models presented in this paper are simple, ide- 
alized models of star-forming cores that are highly pa- 
rameterized. We justified the choice of specific parame- 
terizations and parameter values with theoretical and/or 
observational constraints in most cases. The one notable 
exception is our choice to maximize the momentum en- 
trainment efficiency between the jet/wind and ambient 
medium in i )4.4[ and in this case we discussed the effects 
of varying this parameter. 

In general, we consider our results to be robust to dif- 
ferent values and parameterizations. For example, it is 
the simple presence of a disk in the radiative transfer, 
rather than the exact disk density profile assumed, that 
affects our results since any disk-shaped object will in- 
troduce a similar inclination dependence in calculated 
observational signatures. It is not the exact shape of the 
outflow cavity (assumed here to follow the streamlines of 
the collapse solution and thus be conical at large radii) 
that matters as much as their simple presence in the 
envelope, as outflow cavities of any shape will increase 
Lboi, Thoh E^nd Lfioi/ Lsmm for inclinations that view di- 
rect protostar-t-disk emission through the outflow cavity 
in addition to the far-infrared and millimeter wavelength 
emission from the envelope. It is not the exact choice of 
burst and quiescent accretion rates or the exact condition 
upon which a burst begins that matters as much as it is 
the simple existence of bursts, since a cycle of episodic 
accretion will, in general, shift the models to lower lumi- 
nosities except during bursts and will cause periodic in- 
creases and decreases in evolutionary indicators like Ti,oi 
and Lhoil Lsmm- The details of shape of the model SEDS 
and the comparison to observations will change, but the 
overall results will not. 

To give a quantitative example, motivated by the ex- 
cess of time spent at low-luminosity (~ 0.1 Lq) by model 
5 compared to the observations as seen in Figure [3T] and 
discussed in !j4.51 we constructed an alternate version of 
model 5. Everything remains the same in this alternate 
model except MotoP, the accretion rate from the disk 
onto the protostar, is increased from MjjtoP — Mq 
yr~^ in the quiescent phases between accretion bursts 
to MotoP = 10~^ M0 yr^^, in general agreement with 
Vorobyov & Basu (2005b, 2006). The overaU results of 
this change are minor; the peak in time spent by the 
model at ~ 0.1 Lq decreases by about 25% since the 
nonzero quiescent MjjtoP increases the quiescent phase 
model luminosities, but no substantial changes are in- 
troduced in our conclusions. Figure [32] shows the BLT 
diagram for this alternate version of model 5; while the 
overall distribution of time spent in various bins is simi- 
lar, we note that the location of excluded (white) zones 



changes. These excluded zones should not be consid- 
ered a real effect; slight changes in model parameters 
can move these zones around without changing the over- 
all model conclusions, and increased sampling of the core 
mass function beyond three masses (0.3, 1, and 3 Mq) 
are likely to fill in at least some of them. 




Fig. 32. — Same as Figure l32l except now for the alternate model 

5 with MotoP = 10~^ Mq yr"-"^ in the quiescent accretion phase. 
The grayscale pixels indicate the fraction of total time t he m odel 
spends in each L^ol ~ T^ol bin, calculated from Equation 1101 The 
grayscale is displayed in a logarithmic stretch with the scaling 
chosen to emphasize the full extent of the models in Li,oi ~ T^ol 
space. The mapping between grayscale and fraction of total time 
is indicated in the legend. The class boundaries in Ti,oi are taken 
from Chen et al. (1995). The thick dashed line shows the ZAMS 
(D'Antona & MazziteUi 1994) from 0.1 to 2.0 Mq. The colored 
symbols show the Young Stellar Objects from Evans et al. (2009) 
in this diagram; the colors and symbols hold the same meaning as 
in Figure [3 

A more likely cause of uncertainty than the model pa- 
rameterizations is in the weighting of the different initial 
mass cores by the CMF in order to calculate the frac- 
tions of total time spent in different Lhoj and T^oi bins 
for comparison to observations. The three initial masses 
(0.3, 1, and 3 M©) were chosen both by YE05 and by us 
because they adequately sample both sides of the peak 
of the CMF (~ 1 Mq). However, the exact shape of 
the CMF and thus the relative numbers of different mass 
cores remains a significant unknown, especially at the 
low end where many studies suffer from incompleteness 
effects. The individual model results are robust to differ- 
ent parameterizations, but the combined comparison to 
observations is significnatly more uncertain. This com- 
parison must be revisited as better studies of the CMF 
become available with new instruments such as SCUBA- 
II (Ward-Thompson et al. 2007). 

Finally, we note that the collapse of the core to form a 
protostar follows the inside-out collapse of static, singu- 
lar isothermal spheres first calculated by Shu (1977) and 
extended by TSC84 to include rotation. Other collapse 
solutions that take into account nonzero initial velocities 
(e.g.. Hunter 1977; Fatuzzo et al. 2004) and magnetic 
fields (e.g., Li & Shu 1997) tend to increase the accretion 
rate and would thus worsen the luminosity problem. On 
the other hand, Vorobyov & Basu (2005a) showed that 
a finite mass reservoir will create a phase of terminally 
declining accretion rate, an effect included in their col- 
lapse simulations featuring episodic accretion (Vorobyov 

6 Basu 2005b, 2006). More detailed future models that 
follow the exact evolution of the Vorobyov & Basu sim- 
ulations rather than the simple, idealized models pre- 
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sented here will be needed to fully evaluate the necessity 
and ability of episodic accretion to resolve the luminosity 
problem. 

6. SUMMARY 

We have made five modifications to the YE05 evolu- 
tionary model in an effort to bring the model in better 
agreement with observations: (1) we modified the dust 
opacities to include isotropic scattering off dust grains 
fi j4.ip . (2) we included a circumstellar disk directly in 
the radiative transfer f H4.2p . (3) we included a rotation- 
ally flattened envelope density structure following the 
TSC84 solution for the collapse of slowly rotating cores 
(i i4.3p . (4) we inclu ded the effects of mass-loss and out- 
flow cavities (SHIi and (5) we included a simple treat- 
ment of episodic mass accretion based on the simulations 
by Vorobyov & Basu (2005b, 2006; U^ . 

We find that the first four models all affect the model 
predictions but are unable to resolve the long-standing 
luminosity problem. Including a cycle of episodic ac- 
cretion, however, can resolve this problem and bring the 
model predictions in better agreement with observations. 
We find that standard assumptions about the interplay 
between mass accretion and mass loss in our model give 
star formation efficiencies consistent with recent obser- 
vations that compare the core mass function (CMF) to 



the stellar initial mass function (IMF), and that the com- 
bination of episodic accretion and increased inclination 
dependence introduced by the presence of outfiow cavi- 
ties both work to reduce the connection between physical 
Stage and observational Class as calculated by common 
evolutionary indicators. We have outlined future stud- 
ies needed on both observational and modeling fronts in 
order to test the conclusions of this paper that episodic 
accretion is both necessary and sufficient to resolve the 
luminosity problem. 
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TABLE 1 
Fraction of Time in BLT Bins 



Bin 

10 < Tboi < 100 K 
100 < Tioi < 1000 K 
L0, 1000 < Thoi < 10000 K 
L©, Ttoi > 10000 K 

< 10 Lq, 10 < Tjoi < 100 K 

< 10 L0, 100 < Ttoi < 1000 K 

< 10 L0, 1000 < Ttoi < 10000 K 

< 10 Lq, Ttoi > 10000 K 

< 1000 L0, 10 < Thoi < 100 K 



Lboi < 0.1 L0 
Lbol < 0.1 Lq 
Ltoi < 0.1 
Lboi < 0.1 
0.1 < Lb,i 
0.1 < Lboi 
0.1 < Lboi 
0.1 < Lboi 

10 < Lboi ^, _ „„. - 

10 < Lboi < 1000 Lq, 100 < Tboi < 1000 K 
10 < Lboi < 1000 Lq, 1000 < Tboi < 10000 K 
10 < Lboi < 1000 Lq, Tboi > 10000 K 
Lboi > 1000 Lq, 10 < Tboi < 100 K 
Lboi > 1000 Lq, 100 < Tboi < 1000 K 
Lboi > 1000 Lq, 1000 < Tboi < 10000 K 
Lboi > 1000 Lq, Tboi > 10000 K 



)bscrvations 


YE05 Model 


Model 1 


Model 2 


Model 3 


Model 4 


Model 5 


0% 


5.4% 


6.0% 


6.0% 


5.5% 


9.4% 


19.6% 


5.4% 


0% 


0% 


0% 


0% 


0.2% 


5.5% 


1.8% 


0% 


0% 


0% 


0% 


0.6% 


0.7% 


0% 


0% 


0% 


0% 


0% 


0% 


0% 


32.1% 


11.5% 


12.4% 


20.3% 


21.3% 


34.4% 


28.4% 


43.8% 


5.8% 


0.9% 


0.6% 


0.6% 


11.8% 


8.5% 


0.9% 


0.1% 


0.1% 


0.8% 


0% 


18.3% 


27.4% 


0% 


0% 


0% 


0% 


0% 


0% 


0% 


2.7% 


1.0% 


25.5% 


26.2% 


18.7% 


2.8% 


1.8% 


11.6% 


58.3% 


51.2% 


41.0% 


53.9% 


1.2% 


0.6% 


1.8% 


17.9% 


3.9% 


5.1% 


0% 


21.2% 


7.1% 


0% 


0% 


0% 


0% 


0% 


0% 
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0% 


0% 


0% 


0% 
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0% 


0% 


0% 


0% 
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0% 


0.1% 
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0% 


0% 


0.2% 



TABLE 2 

Summary of Additions to the YE05 Model 



Model 




Scattering 


2-D Disk 


2-D Envelope 
Density Profile 


Mass-Loss and 
Outfiow Cavities 


Episodic 
Accretion 


Model 1 (i 


4.1 


Y 


N 


N 


N 


N 


Model 2 (1 


4.2 


Y 


Y 


N 


N 


N 


Model 3 11 


4.3 


Y 


Y 


Y 


N 


N 


Model 4 (i 


4.4 


Y 


Y 


Y 


Y 


N 


Model 5 (i 


4.5 


Y 


Y 


Y 


Y 


Y 



TABLE 3 

Fraction of Time in Class/Stage Combinations 





Stage 


Stage I 




Mint/Mtot < 0.5 


Mi„t/Mtot > 0.5 


Tboi Class (Tboi < 70 K) 


31% 


6% 


Class I (70 < Tboi < 650 K) 


12% 


9% 


Class II [Tboi > 650 K) 


10% 


31% 


Lbol/Lsmm Class (Lbol/Lsrnm < 175 K) 


43% 


10% 


Class I (Lboi/Lsmm > 175 K) 


11% 


36% 



